diff --git a/Project.toml b/Project.toml
index 40cb98a..b9c52d1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -13,7 +13,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-MultivariatePolynomials = "0.5"
+MultivariatePolynomials = "0.5.3"
 MutableArithmetics = "1"
 Reexport = "1"
 julia = "1"
diff --git a/src/DynamicPolynomials.jl b/src/DynamicPolynomials.jl
index 5df9ab0..3cc41dd 100644
--- a/src/DynamicPolynomials.jl
+++ b/src/DynamicPolynomials.jl
@@ -72,11 +72,11 @@ function MP.similar_variable(
 ) where {V,M,S}
     return MP.similar_variable(P, S)
 end
-function MP.similar_variable(
-    ::Union{PolyType{V,M},Type{<:PolyType{V,M}}},
-    s::Symbol,
-) where {V,M}
-    return Variable(string(s), V, M)
+function MP.similar_variable(p::PolyType{V,M}, s::Symbol) where {V,M}
+    return Variable(string(s), V, M, isreal(p) ? REAL : COMPLEX)
+end
+function MP.similar_variable(::Type{<:PolyType{V,M}}, s::Symbol) where {V,M}
+    return Variable(string(s), V, M, REAL) # we cannot infer this from the type
 end
 
 include("promote.jl")
diff --git a/src/comp.jl b/src/comp.jl
index d9c0fbd..62e7edd 100644
--- a/src/comp.jl
+++ b/src/comp.jl
@@ -27,14 +27,19 @@ function (==)(
     x::Variable{<:AnyCommutative{CreationOrder}},
     y::Variable{<:AnyCommutative{CreationOrder}},
 )
-    return x.variable_order.order.id == y.variable_order.order.id
+    return x.variable_order.order.id == y.variable_order.order.id &&
+           x.kind == y.kind
 end
 
 function Base.isless(
     x::Variable{<:AnyCommutative{CreationOrder}},
     y::Variable{<:AnyCommutative{CreationOrder}},
 )
-    return isless(y.variable_order.order.id, x.variable_order.order.id)
+    if x.variable_order.order.id == y.variable_order.order.id
+        return isless(y.kind, x.kind)
+    else
+        return isless(y.variable_order.order.id, x.variable_order.order.id)
+    end
 end
 
 # Comparison of Monomial
diff --git a/src/mono.jl b/src/mono.jl
index b3d97c6..02ef0c9 100644
--- a/src/mono.jl
+++ b/src/mono.jl
@@ -113,3 +113,10 @@ function __add_variables!(
     mono.z[map] = tmp
     return
 end
+
+# for efficiency reasons
+function Base.conj(x::Monomial{V,M}) where {V<:Commutative,M}
+    cv = conj.(x.vars)
+    perm = sortperm(cv, rev = true)
+    return Monomial{V,M}(cv[perm], x.z[perm])
+end
diff --git a/src/monomial_vector.jl b/src/monomial_vector.jl
index 991f4b1..b9d3e51 100644
--- a/src/monomial_vector.jl
+++ b/src/monomial_vector.jl
@@ -88,6 +88,28 @@ end
 function MultivariatePolynomials.maxdegree(x::MonomialVector)
     return isempty(x) ? 0 : maximum(sum.(x.Z))
 end
+# Complex-valued degrees for monomial vectors
+for (fun, call, def, ret) in [
+    (:extdegree_complex, :extrema, (0, 0), :((min(v1, v2), max(v1, v2)))),
+    (:mindegree_complex, :minimum, 0, :(min(v1, v2))),
+    (:maxdegree_complex, :maximum, 0, :(max(v1, v2))),
+]
+    eval(quote
+        function MultivariatePolynomials.$fun(x::MonomialVector)
+            isempty(x) && return $def
+            vars = variables(x)
+            @assert(!any(isrealpart, vars) && !any(isimagpart, vars))
+            grouping = isconj.(vars)
+            v1 = $call(sum(z[grouping]) for z in x.Z)
+            grouping = map(!, grouping)
+            v2 = $call(sum(z[grouping]) for z in x.Z)
+            return $ret
+        end
+    end)
+end
+# faster complex-related functions
+Base.isreal(x::MonomialVector) = all(isreal, x.vars)
+Base.conj(x::MonomialVector) = MonomialVector(conj.(x.vars), x.Z)
 
 MP.variables(m::Union{Monomial,MonomialVector}) = m.vars
 
@@ -118,7 +140,11 @@ end
 # TODO replace by MP function
 function _error_for_negative_degree(deg)
     if deg < 0
-        throw(ArgumentError("The degree should be a nonnegative number but the provided degree `$deg` is negative."))
+        throw(
+            ArgumentError(
+                "The degree should be a nonnegative number but the provided degree `$deg` is negative.",
+            ),
+        )
     end
 end
 
diff --git a/src/subs.jl b/src/subs.jl
index 216ca00..f856522 100644
--- a/src/subs.jl
+++ b/src/subs.jl
@@ -9,23 +9,46 @@ end
 
 function fillmap!(
     vals,
-    vars::Vector{<:Variable{<:NonCommutative}},
+    vars::Vector{<:Variable{C}},
     s::MP.Substitution,
-)
-    for j in eachindex(vars)
-        if vars[j] == s.first
-            vals[j] = s.second
+) where {C}
+    # We may assign a complex or real variable to its value (ordinary substitution).
+    # We follow the following rules:
+    # - If a single substitution rule determines the value of a real variable, just substitute it.
+    # - If a single substitution rule determines the value of a complex variable or its conjugate, substitute the appropriate
+    #   value whereever something related to this variable is found (i.e., the complex variable, the conjugate variable, or
+    #   its real or imaginary part)
+    # - If a single substitution rule determines the value of the real or imaginary part of a complex variable alone, then only
+    #   replace the real or imaginary parts if they occur explicitly. Don't do a partial substitution, i.e., `z` with the rule
+    #   `zᵣ => 1` is left alone and not changed into `1 + im*zᵢ`. Even if both the real and imaginary parts are substituted as
+    #   two individual rules (which we don't know of in this method), `z` will not be replaced.
+    if s.first.kind == REAL
+        for j in eachindex(vars)
+            if vars[j] == s.first
+                vals[j] = s.second
+                C == Commutative && break
+            end
+        end
+    else
+        s.first.kind == COMPLEX || throw(
+            ArgumentError(
+                "Substitution with complex variables requires the ordinary_variable in the substitution specification",
+            ),
+        )
+        for j in eachindex(vars)
+            if vars[j].variable_order.order.id ==
+               s.first.variable_order.order.id
+                if vars[j].kind == COMPLEX
+                    vals[j] = s.second
+                elseif vars[j].kind == CONJ
+                    vals[j] = conj(s.second)
+                elseif vars[j].kind == REAL_PART
+                    vals[j] = real(s.second)
+                else
+                    vals[j] = imag(s.second)
+                end
+            end
         end
-    end
-end
-function fillmap!(
-    vals,
-    vars::Vector{<:Variable{<:Commutative}},
-    s::MP.Substitution,
-)
-    j = findfirst(isequal(s.first), vars)
-    if j !== nothing
-        vals[j] = s.second
     end
 end
 function fillmap!(vals, vars, s::MP.AbstractMultiSubstitution)
diff --git a/src/var.jl b/src/var.jl
index b7eb350..3e4fc66 100644
--- a/src/var.jl
+++ b/src/var.jl
@@ -1,20 +1,34 @@
-export Variable, @polyvar, @ncpolyvar
+export Variable, @polyvar, @ncpolyvar, @complex_polyvar
 
-function polyarrayvar(variable_order, monomial_order, prefix, indices...)
+function polyarrayvar(
+    variable_order,
+    monomial_order,
+    complex_kind,
+    prefix,
+    indices...,
+)
     return map(
         i -> Variable(
             "$(prefix)[$(join(i, ","))]",
             variable_order,
             monomial_order,
+            complex_kind,
         ),
         Iterators.product(indices...),
     )
 end
 
-function buildpolyvar(var, variable_order, monomial_order)
+function buildpolyvar(var, variable_order, monomial_order, complex_kind)
     if isa(var, Symbol)
         var,
-        :($(esc(var)) = $Variable($"$var", $variable_order, $monomial_order))
+        :(
+            $(esc(var)) = $Variable(
+                $"$var",
+                $variable_order,
+                $monomial_order,
+                $complex_kind,
+            )
+        )
     else
         isa(var, Expr) || error("Expected $var to be a variable name")
         Base.Meta.isexpr(var, :ref) ||
@@ -28,6 +42,7 @@ function buildpolyvar(var, variable_order, monomial_order)
             $(esc(varname)) = polyarrayvar(
                 $(variable_order),
                 $(monomial_order),
+                $(complex_kind),
                 $prefix,
                 $(esc.(var.args[2:end])...),
             )
@@ -35,11 +50,12 @@ function buildpolyvar(var, variable_order, monomial_order)
     end
 end
 
-function buildpolyvars(args, variable_order, monomial_order)
+function buildpolyvars(args, variable_order, monomial_order, complex_kind)
     vars = Symbol[]
     exprs = []
     for arg in args
-        var, expr = buildpolyvar(arg, variable_order, monomial_order)
+        var, expr =
+            buildpolyvar(arg, variable_order, monomial_order, complex_kind)
         push!(vars, var)
         push!(exprs, expr)
     end
@@ -47,7 +63,7 @@ function buildpolyvars(args, variable_order, monomial_order)
 end
 
 # Inspired from `JuMP.Containers._extract_kw_args`
-function _extract_kw_args(args, variable_order)
+function _extract_kw_args(args, variable_order, complex_kind)
     positional_args = Any[]
     monomial_order = :($(MP.Graded{MP.LexOrder}))
     for arg in args
@@ -56,6 +72,8 @@ function _extract_kw_args(args, variable_order)
                 variable_order = esc(arg.args[2])
             elseif arg.args[1] == :monomial_order
                 monomial_order = esc(arg.args[2])
+            elseif arg.args[1] == :complex
+                complex_kind = arg.args[2] == :true ? COMPLEX : REAL
             else
                 error("Unrecognized keyword argument `$(arg.args[1])`")
             end
@@ -63,22 +81,33 @@ function _extract_kw_args(args, variable_order)
             push!(positional_args, arg)
         end
     end
-    return positional_args, variable_order, monomial_order
+    return positional_args, variable_order, monomial_order, complex_kind
 end
 
 # Variable vector x returned garanteed to be sorted so that if p is built with x then vars(p) == x
 macro polyvar(args...)
-    pos_args, variable_order, monomial_order =
-        _extract_kw_args(args, :($(Commutative{CreationOrder})))
-    vars, exprs = buildpolyvars(pos_args, variable_order, monomial_order)
+    pos_args, variable_order, monomial_order, complex_kind =
+        _extract_kw_args(args, :($(Commutative{CreationOrder})), REAL)
+    vars, exprs =
+        buildpolyvars(pos_args, variable_order, monomial_order, complex_kind)
     return :($(foldl((x, y) -> :($x; $y), exprs, init = :()));
     $(Expr(:tuple, esc.(vars)...)))
 end
 
 macro ncpolyvar(args...)
-    pos_args, variable_order, monomial_order =
-        _extract_kw_args(args, :($(NonCommutative{CreationOrder})))
-    vars, exprs = buildpolyvars(pos_args, variable_order, monomial_order)
+    pos_args, variable_order, monomial_order, complex_kind =
+        _extract_kw_args(args, :($(NonCommutative{CreationOrder})), REAL)
+    vars, exprs =
+        buildpolyvars(pos_args, variable_order, monomial_order, complex_kind)
+    return :($(foldl((x, y) -> :($x; $y), exprs, init = :()));
+    $(Expr(:tuple, esc.(vars)...)))
+end
+
+macro complex_polyvar(args...)
+    pos_args, variable_order, monomial_order, complex_kind =
+        _extract_kw_args(args, :($(Commutative{CreationOrder})), COMPLEX)
+    vars, exprs =
+        buildpolyvars(pos_args, variable_order, monomial_order, complex_kind)
     return :($(foldl((x, y) -> :($x; $y), exprs, init = :()));
     $(Expr(:tuple, esc.(vars)...)))
 end
@@ -114,19 +143,54 @@ function instantiate(::Type{NonCommutative{O}}) where {O}
     return NonCommutative(instantiate(O))
 end
 
+"""
+    ComplexKind
+
+The type used to identify whether a variable is complex-valued. It has the following instances:
+- `REAL`: the variable is real-valued, [`conj`](@ref) and [`real`](@ref) are identities, [`imag`](@ref) is zero.
+- `COMPLEX`: the variable is a declared complex-valued variable with distinct conjugate, real, and imaginary part.
+  [`ordinary_variable`](@ref) is an identity.
+- `CONJ`: the variable is the conjugate of a declared complex-valued variable. [`ordinary_variable`](@ref) is equal to
+  [`conj`](@ref). It will print with an over-bar.
+- `REAL_PART`: the variable is the real part of a complex-valued variable. It is real-valued itself: [`conj`](@ref) and
+  [`real`](@ref) are identities, [`imag`](@ref) is zero. [`ordinary_variable`](@ref) will give back the original complex-valued
+  _declared_ variable from which the real part was extracted. It will print with an `ᵣ` subscript.
+- `IMAG_PART`: the variable is the imaginary part of a declared complex-valued variable (or the negative imaginary part of the
+  conjugate of a declared complex-valued variable). It is real-valued itself: [`conj`](@ref) and [`real`](@ref) are identities,
+  [`imag`](@ref) is zero. [`ordinary_variable`](@ref) will give back the original complex-valued _declared_ variable from which
+  the imaginary part was extracted. It will print with an `ᵢ` subscript.
+"""
+@enum ComplexKind REAL COMPLEX CONJ REAL_PART IMAG_PART
+
 struct Variable{V,M} <: AbstractVariable
     name::String
+    kind::ComplexKind
     variable_order::V
 
-    function Variable{V,M}(name::AbstractString) where {V<:AbstractVariableOrdering,M<:MP.AbstractMonomialOrdering}
-        return new{V,M}(convert(String, name), instantiate(V))
+    function Variable{V,M}(
+        name::AbstractString,
+        kind::ComplexKind = REAL,
+    ) where {V<:AbstractVariableOrdering,M<:MP.AbstractMonomialOrdering}
+        return new{V,M}(convert(String, name), kind, instantiate(V))
     end
+
     function Variable(
         name::AbstractString,
         ::Type{V},
         ::Type{M},
+        kind::ComplexKind = REAL,
     ) where {V<:AbstractVariableOrdering,M<:MP.AbstractMonomialOrdering}
-        return new{V,M}(convert(String, name), instantiate(V))
+        return new{V,M}(convert(String, name), kind, instantiate(V))
+    end
+
+    function Variable(
+        from::Variable{V,M},
+        kind::ComplexKind,
+    ) where {V<:AbstractVariableOrdering,M<:MP.AbstractMonomialOrdering}
+        (from.kind == REAL && kind == REAL) ||
+            (from.kind != REAL && kind != REAL) ||
+            error("Cannot change the complex type of an existing variable")
+        return new{V,M}(from.name, kind, from.variable_order)
     end
 end
 
@@ -150,6 +214,41 @@ MP.ordering(::Type{Variable{V,M}}) where {V,M} = M
 
 iscomm(::Type{Variable{C}}) where {C} = C
 
+Base.isreal(x::Variable{C}) where {C} = x.kind != COMPLEX && x.kind != CONJ
+MP.isrealpart(x::Variable{C}) where {C} = x.kind == REAL_PART
+MP.isimagpart(x::Variable{C}) where {C} = x.kind == IMAG_PART
+MP.isconj(x::Variable{C}) where {C} = x.kind == CONJ
+function MP.ordinary_variable(x::Variable)
+    return x.kind == REAL || x.kind == COMPLEX ? x : Variable(x, COMPLEX)
+end
+
+function Base.conj(x::Variable)
+    if x.kind == COMPLEX
+        return Variable(x, CONJ)
+    elseif x.kind == CONJ
+        return Variable(x, COMPLEX)
+    else
+        return x
+    end
+end
+
+function Base.real(x::Variable)
+    if x.kind == COMPLEX || x.kind == CONJ
+        return Variable(x, REAL_PART)
+    else
+        return x
+    end
+end
+function Base.imag(x::Variable{V,M}) where {V,M}
+    if x.kind == COMPLEX
+        return Variable(x, IMAG_PART)
+    elseif x.kind == CONJ
+        return _Term{V,M,Int}(-1, Monomial(Variable(x, IMAG_PART)))
+    else
+        return MA.Zero()
+    end
+end
+
 function mergevars_to!(
     vars::Vector{PV},
     varsvec::Vector{Vector{PV}},
diff --git a/test/poly.jl b/test/poly.jl
index 80c55a3..cb0d7b2 100644
--- a/test/poly.jl
+++ b/test/poly.jl
@@ -55,10 +55,16 @@
         p.x[2] == x
 
         @inferred DynamicPolynomials.Polynomial(i -> float(i), [x, x * x])
-        @inferred DynamicPolynomials.Polynomial(i -> float(i), MonomialVector([x * x, x]))
+        @inferred DynamicPolynomials.Polynomial(
+            i -> float(i),
+            MonomialVector([x * x, x]),
+        )
         for p in (
             DynamicPolynomials.Polynomial(i -> float(i), [x, x * x]),
-            DynamicPolynomials.Polynomial(i -> float(i), MonomialVector([x * x, x])),
+            DynamicPolynomials.Polynomial(
+                i -> float(i),
+                MonomialVector([x * x, x]),
+            ),
         )
             @test typeof(p) == PTF
             @test p.a == [1.0, 2.0]
@@ -130,5 +136,16 @@
         )
         @test_throws err t(y => 1)
         @test_throws err p(y => 1)
+
+        @complex_polyvar z
+        pc = z^3 + 2real(z) - 7imag(z^4)
+        @test pc(z => 2 + 3im) == 798 + 9im
+        err = ArgumentError(
+            "Substitution with complex variables requires the ordinary_variable in the substitution specification"
+        )
+        @test_throws err pc(conj(z) => 2 - 3im)
+        @test_throws err pc(real(z) => 2)
+        @test_throws err pc(imag(z) => 3)
+        @test real(pc)(z => 2 + 3im) == 798
     end
 end