diff --git a/docs/polytope.md b/docs/polytope.md
index c4958c1..20ccf1a 100644
--- a/docs/polytope.md
+++ b/docs/polytope.md
@@ -54,7 +54,7 @@ Next, let's create the inner box approximation of the polytope.
 
 ```{code-cell} ipython3
 inner_event = polytope.inner_box_approximation()
-fig.add_traces(inner_event.plot("aquamarine"))
+fig.add_traces(inner_event.plot("#2ca02c"))
 fig.show()
 ```
 
@@ -63,6 +63,6 @@ Let's create the outer box approximation of the polytope.
 ```{code-cell} ipython3
 outer_event = polytope.outer_box_approximation()
 outer_event = outer_event - inner_event
-fig.add_traces(outer_event.plot("aqua"))
+fig.add_traces(outer_event.plot("#d62728"))
 fig.show()
 ```
diff --git a/src/random_events/interval.py b/src/random_events/interval.py
index 842833f..29977f6 100644
--- a/src/random_events/interval.py
+++ b/src/random_events/interval.py
@@ -181,6 +181,9 @@ def contained_integers(self) -> int:
 
         yield from range(rounded_lower, rounded_upper + 1)
 
+    def __deepcopy__(self):
+        return self.__class__(self.lower, self.upper, self.left, self.right)
+
 
 class Interval(sigma_algebra.AbstractCompositeSet):
     simple_sets: SortedSet[SimpleInterval]
diff --git a/src/random_events/product_algebra.py b/src/random_events/product_algebra.py
index df441e1..6098916 100644
--- a/src/random_events/product_algebra.py
+++ b/src/random_events/product_algebra.py
@@ -326,6 +326,9 @@ def fill_missing_variables(self, variables: VariableSet):
             if variable not in self:
                 self[variable] = variable.domain
 
+    def __deepcopy__(self):
+        return self.__class__({variable: assignment.__deepcopy__() for variable, assignment in self.items()})
+
 
 class Event(AbstractCompositeSet):
     """
@@ -441,9 +444,9 @@ def bounding_box(self) -> SimpleEvent:
         for variable in self.all_variables:
             for simple_set in self.simple_sets:
                 if variable not in result:
-                    result[variable] = simple_set[variable]
+                    result[variable] = simple_set[variable].__deepcopy__()
                 else:
-                    result[variable] = result[variable].union_with(simple_set[variable])
+                    result[variable] = result[variable].union_with(simple_set[variable].__deepcopy__())
         return result
 
     def plot(self, color="#636EFA") -> Union[List[go.Scatter], List[go.Mesh3d]]:
diff --git a/src/random_events/set.py b/src/random_events/set.py
index 192e8c1..a6ed90d 100644
--- a/src/random_events/set.py
+++ b/src/random_events/set.py
@@ -78,6 +78,9 @@ def as_composite_set(self) -> AbstractCompositeSet:
         return Set(self)
 
 
+    def __deepcopy__(self):
+        return self.__class__(self.value)
+
 class Set(AbstractCompositeSet):
     simple_sets: SetElementContainer
 
diff --git a/src/random_events/sigma_algebra.py b/src/random_events/sigma_algebra.py
index 86fb721..9461bc6 100644
--- a/src/random_events/sigma_algebra.py
+++ b/src/random_events/sigma_algebra.py
@@ -114,6 +114,10 @@ def as_composite_set(self) -> AbstractCompositeSet:
         """
         raise NotImplementedError
 
+    @abstractmethod
+    def __deepcopy__(self):
+        raise NotImplementedError
+
 
 class AbstractCompositeSet(SubclassJSONSerializer):
     """
@@ -455,6 +459,9 @@ def to_json(self) -> Dict[str, Any]:
     def _from_json(cls, data: Dict[str, Any]) -> Self:
         return cls(*[AbstractSimpleSet.from_json(simple_set) for simple_set in data["simple_sets"]])
 
+    def __deepcopy__(self):
+        return self.__class__(*[ss.__deepcopy__() for ss in self.simple_sets])
+
 
 # Type definitions
 if TYPE_CHECKING:
diff --git a/test/test_product_algebra.py b/test/test_product_algebra.py
index bf66efb..bec8995 100644
--- a/test/test_product_algebra.py
+++ b/test/test_product_algebra.py
@@ -119,11 +119,23 @@ def test_bounding_box(self):
         event_1 = SimpleEvent({self.x: closed(0, 1), self.y: SimpleInterval(0, 1)}).as_composite_set()
         event_2 = SimpleEvent({self.x: closed(1, 2), self.y: Interval(SimpleInterval(3, 4))}).as_composite_set()
         event = event_1 | event_2
+
+        event_before_bounding_box = event.__deepcopy__()
         bounding_box = event.bounding_box()
         result = SimpleEvent({self.x: closed(0, 2),
                               self.y: SimpleInterval(0, 1).as_composite_set() | SimpleInterval(3, 4).as_composite_set()})
         self.assertEqual(bounding_box, result)
-
+        self.assertEqual(event, event_before_bounding_box)
+
+    def test_complex_event_bounding_box_with_references(self):
+        event1 = SimpleEvent({self.x: closed(0, 1) |  closed(2, 3),
+                             self.y: closed(0, 1) |  closed(2, 3)}).as_composite_set()
+        event2 = SimpleEvent({self.x: closed(1, 2) |  closed(3, 4),
+                             self.y: closed(1, 2) |  closed(3, 4)}).as_composite_set()
+        event = event1 | event2
+        event_before_bb = event.__deepcopy__()
+        bb = event.bounding_box()
+        self.assertEqual(event, event_before_bb)
 
 class NoneTypeObjectInDifferenceTestCase(unittest.TestCase):
     x: Continuous = Continuous("x")