diff --git a/extreme_estimator/estimator/abstract_estimator.py b/extreme_estimator/estimator/abstract_estimator.py
index 72e7670f2feb119fb7e78a33de421fa1181faec4..7f705474ae76fc3d5989302d5f60f51f83eec0a8 100644
--- a/extreme_estimator/estimator/abstract_estimator.py
+++ b/extreme_estimator/estimator/abstract_estimator.py
@@ -7,6 +7,10 @@ class AbstractEstimator(object):
     DURATION = 'Average duration'
     MAE_ERROR = 'Mean Average Error'
 
+    # For each estimator, we shall define:
+    # - The loss
+    # - The optimization method for each part of the process
+
     def __init__(self, dataset: AbstractDataset):
         self.dataset = dataset
         self.additional_information = dict()
diff --git a/extreme_estimator/estimator/full_estimator.py b/extreme_estimator/estimator/full_estimator.py
new file mode 100644
index 0000000000000000000000000000000000000000..33fe652827cef315620179248d0faaf71e62163b
--- /dev/null
+++ b/extreme_estimator/estimator/full_estimator.py
@@ -0,0 +1,53 @@
+from extreme_estimator.R_fit.gev_fit.abstract_margin_model import AbstractMarginModel
+from extreme_estimator.R_fit.max_stable_fit.abstract_max_stable_model import AbstractMaxStableModel
+from extreme_estimator.estimator.abstract_estimator import AbstractEstimator
+from extreme_estimator.estimator.margin_estimator import SmoothMarginEstimator
+from extreme_estimator.estimator.max_stable_estimator import MaxStableEstimator
+from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
+
+
+class AbstractFullEstimator(AbstractEstimator):
+    pass
+
+
+class SmoothMarginalsThenUnitaryMsp(AbstractFullEstimator):
+
+    def __init__(self, dataset: AbstractDataset, margin_model: AbstractMarginModel,
+                 max_stable_model: AbstractMaxStableModel):
+        super().__init__(dataset)
+        # Instantiate the two associated estimators
+        self.margin_estimator = SmoothMarginEstimator(dataset=dataset, margin_model=margin_model)
+        self.max_stable_estimator = MaxStableEstimator(dataset=dataset, max_stable_model=max_stable_model)
+
+    def _fit(self):
+        # Estimate the margin parameters
+        self.margin_estimator.fit()
+        # Compute the maxima_normalized
+        maxima_normalized = self.margin_estimator.margin_model.get_maxima_normalized(maxima=self.dataset.maxima,
+                                                                                     df_gev_params=self.margin_estimator.df_gev_params)
+        # Update maxima normalized field through the dataset object
+        print(maxima_normalized)
+        self.dataset.maxima_normalized = maxima_normalized
+        # Estimate the max stable parameters
+        self.max_stable_estimator.fit()
+
+
+class FullEstimatorInASingleStep(AbstractFullEstimator):
+    pass
+
+
+class FullEstimatorInASingleStepWithSmoothMarginals(AbstractFullEstimator):
+    """The method of Gaume, check if its method is in a single step or not"""
+    pass
+
+
+class PointwiseAndThenUnitaryMsp(AbstractFullEstimator):
+    pass
+
+
+class StochasticExpectationMaximization(AbstractFullEstimator):
+    pass
+
+
+class INLAgoesExtremes(AbstractFullEstimator):
+    pass
diff --git a/extreme_estimator/estimator/full_msp_estimator.py b/extreme_estimator/estimator/full_msp_estimator.py
deleted file mode 100644
index d62f9acd52aff5f0b7d5813c0a2dfb60fe46494c..0000000000000000000000000000000000000000
--- a/extreme_estimator/estimator/full_msp_estimator.py
+++ /dev/null
@@ -1,17 +0,0 @@
-from extreme_estimator.estimator.abstract_estimator import AbstractEstimator
-
-
-class FullEstimatorInASingleStep(AbstractEstimator):
-    pass
-
-
-class PointwiseAndThenUnitaryMsp(AbstractEstimator):
-    pass
-
-
-class SmoothMarginalsThenUnitaryMsp(AbstractEstimator):
-    pass
-
-
-class StochasticExpectationMaximization(AbstractEstimator):
-    pass
diff --git a/extreme_estimator/estimator/margin_estimator.py b/extreme_estimator/estimator/margin_estimator.py
index 00aaad385385a8d3df58417283e48f8a0b4b69e1..2f78390baeeb8796b1ce697c145884a5e55419f4 100644
--- a/extreme_estimator/estimator/margin_estimator.py
+++ b/extreme_estimator/estimator/margin_estimator.py
@@ -1,10 +1,27 @@
+from extreme_estimator.R_fit.gev_fit.abstract_margin_model import AbstractMarginModel
 from extreme_estimator.estimator.abstract_estimator import AbstractEstimator
+from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
 
 
-class PointWiseMarginEstimator(AbstractEstimator):
-    pass
+class AbstractMarginEstimator(AbstractEstimator):
+
+    def __init__(self, dataset: AbstractDataset):
+        super().__init__(dataset)
+        assert dataset.temporal_observations.df_maxima is not None
 
 
-class SmoothMarginEstimator(AbstractEstimator):
-    # with different type of marginals: cosntant, linear....
+class PointWiseMarginEstimator(AbstractMarginEstimator):
     pass
+
+
+class SmoothMarginEstimator(AbstractMarginEstimator):
+    """# with different type of marginals: cosntant, linear...."""
+
+    def __init__(self, dataset: AbstractDataset, margin_model: AbstractMarginModel):
+        super().__init__(dataset)
+        self.margin_model = margin_model
+        self.df_gev_params = None
+
+    def _fit(self):
+        self.df_gev_params = self.margin_model.fitmargin(maxima=self.dataset.maxima,
+                                                         coord=self.dataset.coord)
diff --git a/extreme_estimator/estimator/unitary_msp_estimator.py b/extreme_estimator/estimator/max_stable_estimator.py
similarity index 78%
rename from extreme_estimator/estimator/unitary_msp_estimator.py
rename to extreme_estimator/estimator/max_stable_estimator.py
index b66982efbf61b851f8f791f901c27ae09bbfba67..22d1fbaeca4a2b8916b9e804724af88af664f37f 100644
--- a/extreme_estimator/estimator/unitary_msp_estimator.py
+++ b/extreme_estimator/estimator/max_stable_estimator.py
@@ -4,7 +4,7 @@ from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
 import numpy as np
 
 
-class MaxStableEstimator(AbstractEstimator):
+class AbstractMaxStableEstimator(AbstractEstimator):
 
     def __init__(self, dataset: AbstractDataset, max_stable_model: AbstractMaxStableModel):
         super().__init__(dataset=dataset)
@@ -12,9 +12,14 @@ class MaxStableEstimator(AbstractEstimator):
         # Fit parameters
         self.max_stable_params_fitted = None
 
+
+class MaxStableEstimator(AbstractMaxStableEstimator):
+
     def _fit(self):
-        self.max_stable_params_fitted = self.max_stable_model.fitmaxstab(maxima=self.dataset.maxima,
-                                                                         coord=self.dataset.coord)
+        assert self.dataset.maxima_normalized is not None
+        self.max_stable_params_fitted = self.max_stable_model.fitmaxstab(
+            maxima_normalized=self.dataset.maxima_normalized,
+            coord=self.dataset.coord)
 
     def _error(self, true_max_stable_params: dict):
         absolute_errors = {param_name: np.abs(param_true_value - self.max_stable_params_fitted[param_name])
diff --git a/test/__init__.py b/test/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/test/extreme_estimator/__init__.py b/test/extreme_estimator/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/test/extreme_estimator/test_full_estimators.py b/test/extreme_estimator/test_full_estimators.py
new file mode 100644
index 0000000000000000000000000000000000000000..075e56bba37df0a1fe1b193a9eb0aa014a7a3096
--- /dev/null
+++ b/test/extreme_estimator/test_full_estimators.py
@@ -0,0 +1,40 @@
+import unittest
+
+from extreme_estimator.estimator.full_estimator import FullEstimatorInASingleStep, \
+    FullEstimatorInASingleStepWithSmoothMarginals, SmoothMarginalsThenUnitaryMsp
+from spatio_temporal_dataset.dataset.simulation_dataset import MarginDataset
+from spatio_temporal_dataset.spatial_coordinates.generated_coordinates import CircleCoordinatesRadius1
+from test.extreme_estimator.test_margin_estimators import TestMarginEstimators
+from test.extreme_estimator.test_max_stable_estimators import TestMaxStableEstimators
+from itertools import product
+
+
+class TestFullEstimators(unittest.TestCase):
+    DISPLAY = False
+    FULL_ESTIMATORS = [SmoothMarginalsThenUnitaryMsp]
+
+    def setUp(self):
+        super().setUp()
+        self.spatial_coord = CircleCoordinatesRadius1.from_nb_points(nb_points=5, max_radius=1)
+        self.max_stable_models = TestMaxStableEstimators.load_max_stable_models()
+        self.margin_models = TestMarginEstimators.load_margin_models()
+
+    def test_full_estimators(self):
+        print(self.margin_models, self.max_stable_models)
+        for margin_model, max_stable_model in product(self.margin_models, self.max_stable_models):
+            dataset = MarginDataset.from_sampling(nb_obs=10, margin_model=margin_model,
+                                                  spatial_coordinates=self.spatial_coord)
+
+            for estimator_class in self.FULL_ESTIMATORS:
+                estimator = estimator_class(dataset=dataset, margin_model=margin_model,
+                                            max_stable_model=max_stable_model)
+                estimator.fit()
+                if self.DISPLAY:
+                    print(type(margin_model))
+                    print(dataset.df_dataset.head())
+                    print(estimator.additional_information)
+            self.assertTrue(True)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/test/extreme_estimator/test_gev_fit.py b/test/extreme_estimator/test_gev_fit.py
deleted file mode 100644
index 4f9fd23e1931afe1a96d7f3b8b5c647796613c58..0000000000000000000000000000000000000000
--- a/test/extreme_estimator/test_gev_fit.py
+++ /dev/null
@@ -1,37 +0,0 @@
-import unittest
-import numpy as np
-
-from extreme_estimator.R_fit.gev_fit.gev_mle_fit import GevMleFit
-from extreme_estimator.R_fit.utils import get_loaded_r
-
-
-class TestGevFit(unittest.TestCase):
-
-    def test_unitary_mle_gev_fit(self):
-        r = get_loaded_r()
-        r("""
-        set.seed(42)
-        N <- 50
-        loc = 0; scale = 1; shape <- 1
-        x_gev <- rgev(N, loc = loc, scale = scale, shape = shape)
-        start_loc = 0; start_scale = 1; start_shape = 1
-        """)
-        # Get the MLE estimator
-        estimator = GevMleFit(x_gev=np.array(r['x_gev']),
-                              start_loc=np.float(r['start_loc'][0]),
-                              start_scale=np.float(r['start_scale'][0]),
-                              start_shape=np.float(r['start_shape'][0]))
-        # Compare the MLE estimated parameters to the reference
-        mle_params_estimated = estimator.mle_params
-        mle_params_ref = {'loc': 0.0219, 'scale': 1.0347, 'shape': 0.8290}
-        for key in mle_params_ref.keys():
-            self.assertAlmostEqual(mle_params_ref[key], mle_params_estimated[key], places=3)
-
-    # def test_unitary_frechet_unitary_transformation(self):
-    #     data = np.array([0.0, 1.0])
-    #     # frechet_unitary_transformation()
-    #     self.assertTrue(False)
-
-
-if __name__ == '__main__':
-    unittest.main()
\ No newline at end of file
diff --git a/test/extreme_estimator/test_margin_estimators.py b/test/extreme_estimator/test_margin_estimators.py
new file mode 100644
index 0000000000000000000000000000000000000000..6ea4da1dd564c08caed183241550730c6ff566b7
--- /dev/null
+++ b/test/extreme_estimator/test_margin_estimators.py
@@ -0,0 +1,63 @@
+import unittest
+
+import numpy as np
+
+from extreme_estimator.R_fit.gev_fit.abstract_margin_model import ConstantMarginModel
+from extreme_estimator.R_fit.gev_fit.gev_mle_fit import GevMleFit
+from extreme_estimator.R_fit.utils import get_loaded_r
+from extreme_estimator.estimator.margin_estimator import SmoothMarginEstimator
+from spatio_temporal_dataset.dataset.simulation_dataset import MarginDataset
+from spatio_temporal_dataset.spatial_coordinates.generated_coordinates import CircleCoordinatesRadius1
+
+
+class TestMarginEstimators(unittest.TestCase):
+    DISPLAY = False
+    MARGIN_TYPES = [ConstantMarginModel]
+    MARGIN_ESTIMATORS = [SmoothMarginEstimator]
+
+    def test_unitary_mle_gev_fit(self):
+        r = get_loaded_r()
+        r("""
+        set.seed(42)
+        N <- 50
+        loc = 0; scale = 1; shape <- 1
+        x_gev <- rgev(N, loc = loc, scale = scale, shape = shape)
+        start_loc = 0; start_scale = 1; start_shape = 1
+        """)
+        # Get the MLE estimator
+        estimator = GevMleFit(x_gev=np.array(r['x_gev']),
+                              start_loc=np.float(r['start_loc'][0]),
+                              start_scale=np.float(r['start_scale'][0]),
+                              start_shape=np.float(r['start_shape'][0]))
+        # Compare the MLE estimated parameters to the reference
+        mle_params_estimated = estimator.mle_params
+        mle_params_ref = {'loc': 0.0219, 'scale': 1.0347, 'shape': 0.8290}
+        for key in mle_params_ref.keys():
+            self.assertAlmostEqual(mle_params_ref[key], mle_params_estimated[key], places=3)
+
+    def setUp(self):
+        super().setUp()
+        self.spatial_coord = CircleCoordinatesRadius1.from_nb_points(nb_points=5, max_radius=1)
+        self.margin_models = self.load_margin_models()
+
+    @classmethod
+    def load_margin_models(cls):
+        return [margin_class() for margin_class in cls.MARGIN_TYPES]
+
+    def test_dependency_estimators(self):
+        for margin_model in self.margin_models:
+            dataset = MarginDataset.from_sampling(nb_obs=10, margin_model=margin_model,
+                                                  spatial_coordinates=self.spatial_coord)
+
+            for estimator_class in self.MARGIN_ESTIMATORS:
+                estimator = estimator_class(dataset=dataset, margin_model=margin_model)
+                estimator.fit()
+                if self.DISPLAY:
+                    print(type(margin_model))
+                    print(dataset.df_dataset.head())
+                    print(estimator.additional_information)
+            self.assertTrue(True)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/test/extreme_estimator/test_max_stable_estimators.py b/test/extreme_estimator/test_max_stable_estimators.py
new file mode 100644
index 0000000000000000000000000000000000000000..f8365192bdee511618a1a8a42f39bd07e7203e53
--- /dev/null
+++ b/test/extreme_estimator/test_max_stable_estimators.py
@@ -0,0 +1,51 @@
+import unittest
+
+from extreme_estimator.R_fit.max_stable_fit.abstract_max_stable_model import \
+    AbstractMaxStableModelWithCovarianceFunction, CovarianceFunction
+from extreme_estimator.R_fit.max_stable_fit.max_stable_models import Smith, BrownResnick, Schlather, \
+    Geometric, ExtremalT, ISchlather
+from extreme_estimator.estimator.max_stable_estimator import MaxStableEstimator
+from spatio_temporal_dataset.dataset.simulation_dataset import MaxStableDataset
+from spatio_temporal_dataset.spatial_coordinates.generated_coordinates import CircleCoordinatesRadius1
+
+
+class TestMaxStableEstimators(unittest.TestCase):
+    DISPLAY = False
+    MAX_STABLE_TYPES = [Smith, BrownResnick, Schlather, Geometric, ExtremalT, ISchlather]
+    MAX_STABLE_ESTIMATORS = [MaxStableEstimator]
+
+    def setUp(self):
+        super().setUp()
+        self.spatial_coord = CircleCoordinatesRadius1.from_nb_points(nb_points=5, max_radius=1)
+        self.max_stable_models = self.load_max_stable_models()
+
+    @classmethod
+    def load_max_stable_models(cls):
+        # Load all max stable model
+        max_stable_models = []
+        for max_stable_class in cls.MAX_STABLE_TYPES:
+            if issubclass(max_stable_class, AbstractMaxStableModelWithCovarianceFunction):
+                max_stable_models.extend([max_stable_class(covariance_function=covariance_function)
+                                               for covariance_function in CovarianceFunction])
+            else:
+                max_stable_models.append(max_stable_class())
+        return max_stable_models
+
+    def test_max_stable_estimators(self):
+        for max_stable_model in self.max_stable_models:
+            dataset = MaxStableDataset.from_sampling(nb_obs=10,
+                                                     max_stable_model=max_stable_model,
+                                                     spatial_coordinates=self.spatial_coord)
+
+            for estimator_class in self.MAX_STABLE_ESTIMATORS:
+                estimator = estimator_class(dataset=dataset, max_stable_model=max_stable_model)
+                estimator.fit()
+                if self.DISPLAY:
+                    print(type(max_stable_model))
+                    print(dataset.df_dataset.head())
+                    print(estimator.additional_information)
+        self.assertTrue(True)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/test/extreme_estimator/test_max_stable_fit.py b/test/extreme_estimator/test_max_stable_fit.py
deleted file mode 100644
index 9ea78042661010db918ccaec004a51e6840732eb..0000000000000000000000000000000000000000
--- a/test/extreme_estimator/test_max_stable_fit.py
+++ /dev/null
@@ -1,38 +0,0 @@
-import unittest
-
-from extreme_estimator.R_fit.max_stable_fit.abstract_max_stable_model import \
-    AbstractMaxStableModelWithCovarianceFunction, CovarianceFunction
-from extreme_estimator.R_fit.max_stable_fit.max_stable_models import Smith, BrownResnick, Schlather, \
-    Geometric, ExtremalT, ISchlather
-from spatio_temporal_dataset.dataset.simulation_dataset import SimulatedDataset
-from spatio_temporal_dataset.spatial_coordinates.generated_coordinates import CircleCoordinatesRadius1
-
-
-class TestMaxStableFit(unittest.TestCase):
-    MAX_STABLE_CLASSES = [Smith, BrownResnick, Schlather, Geometric, ExtremalT, ISchlather]
-
-    def setUp(self):
-        super().setUp()
-        self.spatial_coord = CircleCoordinatesRadius1.from_nb_points(nb_points=5, max_radius=1)
-        self.max_stable_models = []
-        for max_stable_class in self.MAX_STABLE_CLASSES:
-            if issubclass(max_stable_class, AbstractMaxStableModelWithCovarianceFunction):
-                self.max_stable_models.extend([max_stable_class(covariance_function=covariance_function)
-                                               for covariance_function in CovarianceFunction])
-            else:
-                self.max_stable_models.append(max_stable_class())
-
-    def test_sampling_fit_with_models(self, display=False):
-        for max_stable_model in self.max_stable_models:
-            dataset = SimulatedDataset.from_max_stable_sampling(nb_obs=10, max_stable_model=max_stable_model,
-                                                                spatial_coordinates=self.spatial_coord)
-            fitted_values = max_stable_model.fitmaxstab(maxima=dataset.maxima, coord=dataset.coord)
-            if display:
-                print(type(max_stable_model))
-                print(dataset.df_dataset.head())
-                print(fitted_values)
-        self.assertTrue(True)
-
-
-if __name__ == '__main__':
-    unittest.main()
diff --git a/test/extreme_estimator/test_pipeline.py b/test/extreme_estimator/test_pipeline.py
deleted file mode 100644
index d3357fefce861bc131732b076823b06bcdd53ada..0000000000000000000000000000000000000000
--- a/test/extreme_estimator/test_pipeline.py
+++ /dev/null
@@ -1,52 +0,0 @@
-import unittest
-import pandas as pd
-
-class TestPipeline(unittest.TestCase):
-
-    def main_pipeline(self):
-        # Select a type of marginals (either spatial, spatio temporal, temporal)
-        #  this will define the dimension of the climatic space of interest
-        pass
-        # Select the max stable
-
-        #  Define an optimization process
-        # The algo: In 1 time, in 2 times, ..., or more complex patterns
-        # This algo have at least main procedures (that might be repeated several times)
-
-        # For each procedure, we shall define:
-        # - The loss
-        # - The optimization method for each part of the process
-
-    def blanchet_smooth_pipeline(self):
-        pass
-        # Spatial marginal
-
-        # NO MAX STABLE
-
-        # Procedure:
-        # Optimization of a single likelihood process that sums up the likelihood of all the terms.
-
-    def padoan_extreme_pipeline(self):
-        pass
-        # Spatial marginal
-
-        # todo: question, when we are optimizing the full Pairwise loss, are we just optimization the relations ?
-        # or ideally do we need to add the term of order 1
-
-    def gaume(self):
-        # Combining the 2
-        pass
-
-
-
-
-    def test_pipeline_spatial(self):
-        pass
-        # Sample from a
-
-        # Fit the spatio temporal experiment margin
-
-        # Fit the max stable process
-
-if __name__ == '__main__':
-    unittest.main()
\ No newline at end of file
diff --git a/test/spatio_temporal_dataset/__init__.py b/test/spatio_temporal_dataset/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/test/spatio_temporal_dataset/test_temporal_observations.py b/test/spatio_temporal_dataset/test_temporal_observations.py
new file mode 100644
index 0000000000000000000000000000000000000000..0ffd73059dea5780e4d9e824fe67a206d8fa1988
--- /dev/null
+++ b/test/spatio_temporal_dataset/test_temporal_observations.py
@@ -0,0 +1,15 @@
+import unittest
+
+from spatio_temporal_dataset.spatial_coordinates.alps_station_2D_coordinates import \
+    AlpsStation2DCoordinatesBetweenZeroAndOne
+from spatio_temporal_dataset.spatial_coordinates.alps_station_3D_coordinates import \
+    AlpsStation3DCoordinatesWithAnisotropy
+from spatio_temporal_dataset.spatial_coordinates.generated_coordinates import CircleCoordinatesRadius1
+
+
+class TestTemporalObservations(unittest.TestCase):
+
+    DISPLAY = False
+
+if __name__ == '__main__':
+    unittest.main()