diff --git a/extreme_fit/function/param_function/param_function.py b/extreme_fit/function/param_function/param_function.py
index ceb79dec8518f9dc64363dd7263889abeca16715..f59c925d6a268876b672b2e6b4f32a44980d975d 100644
--- a/extreme_fit/function/param_function/param_function.py
+++ b/extreme_fit/function/param_function/param_function.py
@@ -59,8 +59,11 @@ class PolynomialParamFunction(AbstractParamFunction):
 
     def get_param_value(self, coordinate: np.ndarray) -> float:
         gev_param_value = 0
-        for dim, max_degree in self.dim_and_degree:
-            for degree in range(max_degree+1):
+        for i, (dim, max_degree) in enumerate(self.dim_and_degree):
+            # Add intercept only once
+            add_intercept = i == 0
+            first_degree = 0 if add_intercept else 1
+            for degree in range(first_degree, max_degree+1):
                 polynomial_coef = self.coef.dim_to_polynomial_coef[dim]  # type: PolynomialCoef
                 polynomial_coef_value = polynomial_coef.idx_to_coef[degree]
                 gev_param_value += polynomial_coef_value * np.power(coordinate[dim], degree)
diff --git a/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py b/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py
index d2f301cb5a53a1954fc6d069c1efb912eedeeaf5..2770f43a75f378dec3a57e4acde779bd922a109e 100644
--- a/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py
+++ b/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py
@@ -40,9 +40,9 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
             assert len(data) == len(df_coordinates_temp.values), 'len(data)={} != len(temp)={}'.format(len(data),
                                                                                                        len(
                                                                                                            df_coordinates_temp.values))
-            x = ro.FloatVector(data)
         else:
-            x = ro.Matrix(data)
+            data = data.flatten()
+        x = ro.FloatVector(data)
         if self.params_class is GevParams:
             if self.fit_method == MarginFitMethod.is_mev_gev_fit:
                 return self.ismev_gev_fit(x, df_coordinates_temp)
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
index 124093a4cf8bf85f20b513adc6fd770a7f850a6f..607f5fc14a9f87f0f99d7638f9c3af3311da31fb 100644
--- a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_mle_extremes.py
@@ -20,7 +20,6 @@ class ResultFromMleExtremes(AbstractResultFromExtremes):
     @property
     def margin_coef_ordered_dict(self):
         values = self.name_to_value['results']
-        print(values)
         d = self.get_python_dictionary(values)
         if 'par' in d:
             values = {i: param for i, param in enumerate(np.array(d['par']))}
diff --git a/extreme_fit/model/result_from_model_fit/utils.py b/extreme_fit/model/result_from_model_fit/utils.py
index d6a7367c58f012b6d3683c696ca4bb3c926e19ae..844639a5ece8bd899ce61fddd88a90e413ec7132 100644
--- a/extreme_fit/model/result_from_model_fit/utils.py
+++ b/extreme_fit/model/result_from_model_fit/utils.py
@@ -33,7 +33,9 @@ def get_margin_coef_ordered_dict(param_name_to_dims, mle_values, type_for_mle="G
                 coef_dict[coef_name] = mle_values[i]
                 i += 1
             else:
-                for dim, max_degree in dims:
+                # We assume that time was the first parameter
+                inverted_dims = dims[::-1]
+                for dim, max_degree in inverted_dims:
                     coordinate_name = dim_to_coordinate_name[dim]
                     coef_template = LinearCoef.coef_template_str(param_name, coordinate_name)
                     for j in range(1, max_degree + 1):
diff --git a/test/test_extreme_fit/test_estimator/test_gev_spatio_temporal_extremes_mle.py b/test/test_extreme_fit/test_estimator/test_gev_spatio_temporal_extremes_mle.py
index c28f3503c598410b76a2d0b0769dfeb6c5f7fbf2..15553a9428b93934ad06bae38d53e0401ed5f60b 100644
--- a/test/test_extreme_fit/test_estimator/test_gev_spatio_temporal_extremes_mle.py
+++ b/test/test_extreme_fit/test_estimator/test_gev_spatio_temporal_extremes_mle.py
@@ -21,6 +21,7 @@ from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoo
 from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_coordinates import \
     AbstractTemporalCoordinates
 from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
+from spatio_temporal_dataset.slicer.split import Split
 from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_temporal_observations import \
     AbstractSpatioTemporalObservations
 from test.test_projects.test_contrasting.test_two_fold_fit import load_two_fold_fit
@@ -43,26 +44,27 @@ class TestGevTemporalQuadraticExtremesMle(unittest.TestCase):
         estimator = model_fit.estimator_fold_1
         return estimator
 
-    # def test_location_spatio_temporal_linearity(self):
-    #     # Create estimator
-    #     # estimator = fitted_linear_margin_estimator(model_class,
-    #     #                                            self.coordinates, self.dataset,
-    #     #                                            starting_year=0,
-    #     #                                            fit_method=self.fit_method)
-    #     # # Checks that parameters returned are indeed different
-    #     # mle_params_estimated_year1 = estimator.function_from_fit.get_params(np.array([1])).to_dict()
-    #     # mle_params_estimated_year3 = estimator.function_from_fit.get_params(np.array([21])).to_dict()
-    #     # mle_params_estimated_year5 = estimator.function_from_fit.get_params(np.array([41])).to_dict()
-    #     # self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
-    #     # self.assertNotEqual(mle_params_estimated_year3, mle_params_estimated_year5)
-    #     # # Assert the relationship for the location is indeed quadratic
-    #     # diff1 = mle_params_estimated_year1[quadratic_param] - mle_params_estimated_year3[quadratic_param]
-    #     # diff2 = mle_params_estimated_year3[quadratic_param] - mle_params_estimated_year5[quadratic_param]
-    #     # self.assertNotAlmostEqual(diff1, diff2)
-    #     estimator = self.get_estimator_fitted(NonStationaryLocationSpatioTemporalLinearityModel)
-    #     # 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())
+    def test_location_spatio_temporal_linearity(self):
+
+        # Create estimator
+        # estimator = fitted_linear_margin_estimator(model_class,
+        #                                            self.coordinates, self.dataset,
+        #                                            starting_year=0,
+        #                                            fit_method=self.fit_method)
+        # # Checks that parameters returned are indeed different
+        # mle_params_estimated_year1 = estimator.function_from_fit.get_params(np.array([1])).to_dict()
+        # mle_params_estimated_year3 = estimator.function_from_fit.get_params(np.array([21])).to_dict()
+        # mle_params_estimated_year5 = estimator.function_from_fit.get_params(np.array([41])).to_dict()
+        # self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
+        # self.assertNotEqual(mle_params_estimated_year3, mle_params_estimated_year5)
+        # # Assert the relationship for the location is indeed quadratic
+        # diff1 = mle_params_estimated_year1[quadratic_param] - mle_params_estimated_year3[quadratic_param]
+        # diff2 = mle_params_estimated_year3[quadratic_param] - mle_params_estimated_year5[quadratic_param]
+        # self.assertNotAlmostEqual(diff1, diff2)
+        estimator = self.get_estimator_fitted(NonStationaryLocationSpatioTemporalLinearityModel)
+        # Assert that indicators are correctly computed
+        self.assertAlmostEqual(estimator.result_from_model_fit.nllh, estimator.nllh(split=estimator.train_split))
+        # self.assertAlmostEqual(estimator.result_from_model_fit.aic, estimator.aic())
     #     self.assertAlmostEqual(estimator.result_from_model_fit.bic, estimator.bic())