diff --git a/extreme_fit/function/margin_function/polynomial_margin_function.py b/extreme_fit/function/margin_function/polynomial_margin_function.py
index 325ac6335f2ce6cc390a2a2813b22019298782b4..06bbeac5fc69e29483bb43c72635eaf1177fcdcd 100644
--- a/extreme_fit/function/margin_function/polynomial_margin_function.py
+++ b/extreme_fit/function/margin_function/polynomial_margin_function.py
@@ -12,20 +12,20 @@ from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoo
 
 class PolynomialMarginFunction(LinearMarginFunction):
 
-    def __init__(self, coordinates: AbstractCoordinates, param_name_to_dim_and_degree: Dict[str, List[Tuple[int, int]]],
+    def __init__(self, coordinates: AbstractCoordinates, param_name_to_dim_and_max_degree: Dict[str, List[Tuple[int, int]]],
                  param_name_to_coef: Dict[str, PolynomialAllCoef], starting_point: Union[None, int] = None,
                  params_class: type = GevParams):
         param_name_to_dims = {}
-        for param_name in param_name_to_dim_and_degree.keys():
-            dims = [c[0] for c in param_name_to_dim_and_degree[param_name]]
+        for param_name in param_name_to_dim_and_max_degree.keys():
+            dims = [c[0] for c in param_name_to_dim_and_max_degree[param_name]]
             param_name_to_dims[param_name] = dims
-        self.param_name_to_dim_and_degree = param_name_to_dim_and_degree
+        self.param_name_to_dim_and_max_degree = param_name_to_dim_and_max_degree
         super().__init__(coordinates, param_name_to_dims, param_name_to_coef, starting_point, params_class)
 
     COEF_CLASS = PolynomialAllCoef
 
     def load_specific_param_function(self, param_name):
-        return PolynomialParamFunction(dim_and_degree=self.param_name_to_dim_and_degree[param_name],
+        return PolynomialParamFunction(dim_and_degree=self.param_name_to_dim_and_max_degree[param_name],
                                        coef=self.param_name_to_coef[param_name])
 
     def get_params(self, coordinate: np.ndarray, is_transformed: bool = True) -> GevParams:
@@ -34,14 +34,14 @@ class PolynomialMarginFunction(LinearMarginFunction):
     @classmethod
     def from_coef_dict(cls, coordinates: AbstractCoordinates, param_name_to_dims: Dict[str, List[Tuple[int, int]]],
                        coef_dict: Dict[str, float], starting_point: Union[None, int] = None):
-        param_name_to_dim_and_degree = param_name_to_dims
+        param_name_to_dim_and_max_degree = param_name_to_dims
         assert cls.COEF_CLASS is not None, 'a COEF_CLASS class attributes needs to be defined'
         param_name_to_coef = {}
         for param_name in GevParams.PARAM_NAMES:
-            dims = param_name_to_dim_and_degree.get(param_name, [])
+            dims = param_name_to_dim_and_max_degree.get(param_name, [])
             coef = cls.COEF_CLASS.from_coef_dict(coef_dict=coef_dict, param_name=param_name,
                                                  dims=dims,
                                                  coordinates=coordinates)
             param_name_to_coef[param_name] = coef
-        return cls(coordinates, param_name_to_dim_and_degree, param_name_to_coef, starting_point)
+        return cls(coordinates, param_name_to_dim_and_max_degree, param_name_to_coef, starting_point)
 
diff --git a/extreme_fit/function/param_function/polynomial_coef.py b/extreme_fit/function/param_function/polynomial_coef.py
index df8d2e89386d8b2b50defd0ca25e2a49480b78bd..e5784bd10811efd51ee3b451d8b11a118ffc6619 100644
--- a/extreme_fit/function/param_function/polynomial_coef.py
+++ b/extreme_fit/function/param_function/polynomial_coef.py
@@ -20,6 +20,10 @@ class PolynomialCoef(AbstractCoef):
     def compute_default_value(self, idx):
         return self.default_value / idx
 
+    @property
+    def nb_params(self):
+        return self.max_degree + 1
+
 
 class PolynomialAllCoef(LinearCoef):
 
@@ -28,6 +32,13 @@ class PolynomialAllCoef(LinearCoef):
         self.dim_to_polynomial_coef = dim_to_polynomial_coef
         self._intercept = intercept
 
+    @property
+    def nb_params(self):
+        if self.dim_to_polynomial_coef is None:
+            return 1
+        else:
+            return sum([c.nb_params for c in self.dim_to_polynomial_coef.values()])
+
     @property
     def intercept(self) -> float:
         if self._intercept is not None:
diff --git a/extreme_fit/model/margin_model/polynomial_margin_model.py b/extreme_fit/model/margin_model/polynomial_margin_model.py
index 8f40068c8492b8a4b624afc419e93f158ee1e2ab..1f32796148f8f5d19b23c4a27c1a8cd02ccef195 100644
--- a/extreme_fit/model/margin_model/polynomial_margin_model.py
+++ b/extreme_fit/model/margin_model/polynomial_margin_model.py
@@ -20,21 +20,20 @@ class PolynomialMarginModel(AbstractTemporalLinearMarginModel):
                          params_initial_fit_bayesian, type_for_MLE, params_class)
         self.max_degree = max_degree
 
-    # @property
-    # def nb_params(self):
-    #     self.margin_function.param_name_to_coef
-    #     return
+    @property
+    def nb_params(self):
+        return sum([c.nb_params for c in self.margin_function.param_name_to_coef.values()])
 
     @cached_property
     def margin_function(self) -> PolynomialMarginFunction:
         return super().margin_function
 
     def load_margin_function(self, param_name_to_list_dim_and_degree=None):
-        param_name_to_polynomial_all_coef = self.param_name_to_polynomial_all_coef(
-            param_name_and_dim_and_degree_to_coef=self.params_sample)
+        param_name_to_polynomial_all_coef = self.param_name_to_polynomial_all_coef(param_name_to_list_dim_and_degree=param_name_to_list_dim_and_degree,
+                                                                                   param_name_and_dim_and_degree_to_default_coef=self.default_params)
         return PolynomialMarginFunction(coordinates=self.coordinates,
                                         param_name_to_coef=param_name_to_polynomial_all_coef,
-                                        param_name_to_dim_and_degree=param_name_to_list_dim_and_degree,
+                                        param_name_to_dim_and_max_degree=param_name_to_list_dim_and_degree,
                                         starting_point=self.starting_point,
                                         params_class=self.params_class)
 
@@ -48,25 +47,36 @@ class PolynomialMarginModel(AbstractTemporalLinearMarginModel):
                     param_name_and_dim_and_degree_to_coef[(param_name, dim, degree)] = default_slope
         return param_name_and_dim_and_degree_to_coef
 
-    def param_name_to_polynomial_all_coef(self, param_name_and_dim_and_degree_to_coef):
+    def param_name_to_polynomial_all_coef(self, param_name_to_list_dim_and_degree,
+                                          param_name_and_dim_and_degree_to_default_coef):
         param_name_to_polynomial_all_coef = {}
-        param_names = list(set([e[0] for e in param_name_and_dim_and_degree_to_coef.keys()]))
+        param_names = list(set([e[0] for e in param_name_and_dim_and_degree_to_default_coef.keys()]))
         for param_name in param_names:
             dim_to_polynomial_coef = {}
-            for dim in self.coordinates.coordinates_dims:
+            for dim, max_degree in param_name_to_list_dim_and_degree.get(param_name, []):
                 degree_to_coef = {}
-                for (param_name_loop, dim_loop, degree), coef in param_name_and_dim_and_degree_to_coef.items():
+                for (param_name_loop, dim_loop, degree), coef in param_name_and_dim_and_degree_to_default_coef.items():
                     if param_name == param_name_loop and dim == dim_loop:
                         degree_to_coef[degree] = coef
-                dim_to_polynomial_coef[dim] = PolynomialCoef(param_name, degree_to_coef=degree_to_coef)
+                # print(degree_to_coef, param_name)
+                # if len(degree_to_coef) == 0:
+                #     degree_to_coef = {0: param_name_and_dim_and_degree_to_default_coef[(param_name, dim, 0)]}
+                polynomial_coef = PolynomialCoef(param_name, degree_to_coef=degree_to_coef)
+                dim_to_polynomial_coef[dim] = polynomial_coef
+            if len(dim_to_polynomial_coef) == 0:
+                intercept = param_name_and_dim_and_degree_to_default_coef[(param_name, 0, 0)]
+                dim_to_polynomial_coef = None
+            else:
+                intercept = None
             polynomial_all_coef = PolynomialAllCoef(param_name=param_name,
-                                                    dim_to_polynomial_coef=dim_to_polynomial_coef)
+                                                    dim_to_polynomial_coef=dim_to_polynomial_coef,
+                                                    intercept=intercept)
             param_name_to_polynomial_all_coef[param_name] = polynomial_all_coef
         return param_name_to_polynomial_all_coef
 
     @property
     def param_name_to_list_for_result(self):
-        return self.margin_function.param_name_to_dim_and_degree
+        return self.margin_function.param_name_to_dim_and_max_degree
 
 
 class NonStationaryQuadraticLocationModel(PolynomialMarginModel):
diff --git a/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_polynomial_extremes_mle.py b/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_polynomial_extremes_mle.py
index 4c20988b90c9cdc070f33d02034ebe46189b054c..54817ce0fc783377186bad6734e0c9266c0ba1b8 100644
--- a/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_polynomial_extremes_mle.py
+++ b/test/test_extreme_fit/test_distribution/test_gev/test_gev_temporal_polynomial_extremes_mle.py
@@ -57,8 +57,8 @@ class TestGevTemporalQuadraticExtremesMle(unittest.TestCase):
         self.assertNotAlmostEqual(diff1, diff2)
         # Assert that indicators are correctly computed
         self.assertAlmostEqual(estimator.result_from_model_fit.nllh, estimator.nllh())
-        # self.assertAlmostEqual(estimator.result_from_model_fit.aic, estimator.aic())
-        # self.assertAlmostEqual(estimator.result_from_model_fit.bic, estimator.bic())
+        self.assertAlmostEqual(estimator.result_from_model_fit.aic, estimator.aic())
+        self.assertAlmostEqual(estimator.result_from_model_fit.bic, estimator.bic())
 
 
 if __name__ == '__main__':