diff --git a/extreme_data/nasa_data/global_mean_temperature.py b/extreme_data/nasa_data/global_mean_temperature.py
index e8ff3ea1a6d782d2a13d1bde8b11a66fd8ad3be1..8475b218138769ba5b799f2c23c7a3ce0989403d 100644
--- a/extreme_data/nasa_data/global_mean_temperature.py
+++ b/extreme_data/nasa_data/global_mean_temperature.py
@@ -16,5 +16,9 @@ def load_year_to_mean_global_temperature():
     df = pd.read_csv(edf_filepath)
     df = df.astype({'Year': 'float'})
     d = dict(zip(df['Year'], df['Actual Temp']))
+    # Cheat for the last years
+    print('Warning: using fake values for the last few years !!!')
+    for idx in range(1, 5):
+        d[2016 + idx] = d[2016]
     return d
 
diff --git a/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py b/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py
index 476d071e7031262db60c47550f39388c8eebfe56..8a73686eb5e93647c94b6e047762da2074d1e147 100644
--- a/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py
+++ b/extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py
@@ -41,6 +41,7 @@ class LinearMarginEstimator(AbstractMarginEstimator):
     @property
     def df_coordinates_temp(self):
         return self.dataset.coordinates.df_temporal_coordinates_for_fit(split=self.train_split,
+                                                                        temporal_covariate_for_fit=self.margin_model.temporal_covariate_for_fit,
                                                                         starting_point=self.margin_model.starting_point,
                                                                         drop_duplicates=self.margin_model.drop_duplicates)
 
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 2770f43a75f378dec3a57e4acde779bd922a109e..b69a6b4375a2c075c74ea24f58287220fd30bc46 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
@@ -25,13 +25,13 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
                  nb_iterations_for_bayesian_fit=5000,
                  params_initial_fit_bayesian=None,
                  type_for_MLE="GEV",
-                 params_class=GevParams):
-        super().__init__(coordinates, params_user, starting_point, params_class)
+                 params_class=GevParams,
+                 temporal_covariate_for_fit=None):
+        super().__init__(coordinates, params_user, starting_point, params_class, fit_method, temporal_covariate_for_fit)
         self.type_for_mle = type_for_MLE
         self.params_initial_fit_bayesian = params_initial_fit_bayesian
         self.nb_iterations_for_bayesian_fit = nb_iterations_for_bayesian_fit
-        assert isinstance(fit_method, MarginFitMethod), fit_method
-        self.fit_method = fit_method
+        assert isinstance(self.fit_method, MarginFitMethod), self.fit_method
 
     def fitmargin_from_maxima_gev(self, data: np.ndarray, df_coordinates_spat: pd.DataFrame,
                                   df_coordinates_temp: pd.DataFrame) -> AbstractResultFromModelFit:
diff --git a/extreme_fit/model/margin_model/parametric_margin_model.py b/extreme_fit/model/margin_model/parametric_margin_model.py
index 4e4b55315a719a133468728b6c06b1de2f650641..771f8896bcff7f6e87dc4b504053bdad62415c01 100644
--- a/extreme_fit/model/margin_model/parametric_margin_model.py
+++ b/extreme_fit/model/margin_model/parametric_margin_model.py
@@ -20,7 +20,8 @@ class ParametricMarginModel(AbstractMarginModel, ABC):
 
     def __init__(self, coordinates: AbstractCoordinates,
                  params_user=None, starting_point=None, params_class=GevParams,
-                 fit_method=MarginFitMethod.spatial_extremes_mle):
+                 fit_method=MarginFitMethod.spatial_extremes_mle,
+                 temporal_covariate_for_fit=None):
         """
         :param starting_point: starting coordinate for the temporal trend
         """
@@ -28,6 +29,7 @@ class ParametricMarginModel(AbstractMarginModel, ABC):
         self.fit_method = fit_method
         self.starting_point = starting_point
         self.drop_duplicates = True
+        self.temporal_covariate_for_fit = temporal_covariate_for_fit
 
     @cached_property
     def margin_function(self) -> ParametricMarginFunction:
diff --git a/extreme_fit/model/margin_model/polynomial_margin_model/polynomial_margin_model.py b/extreme_fit/model/margin_model/polynomial_margin_model/polynomial_margin_model.py
index 0f044d8ba9cbdd661a47318f9a11603ef53acecc..7357cf84684eefa2d32828fd682e162f2bf3b68f 100644
--- a/extreme_fit/model/margin_model/polynomial_margin_model/polynomial_margin_model.py
+++ b/extreme_fit/model/margin_model/polynomial_margin_model/polynomial_margin_model.py
@@ -17,9 +17,10 @@ class PolynomialMarginModel(AbstractTemporalLinearMarginModel):
 
     def __init__(self, coordinates: AbstractCoordinates, params_user=None, starting_point=None,
                  fit_method=MarginFitMethod.extremes_fevd_mle, nb_iterations_for_bayesian_fit=5000,
-                 params_initial_fit_bayesian=None, type_for_MLE="GEV", params_class=GevParams, max_degree=2):
+                 params_initial_fit_bayesian=None, type_for_MLE="GEV", params_class=GevParams, max_degree=2,
+                 temporal_covariate_for_fit=None):
         super().__init__(coordinates, params_user, starting_point, fit_method, nb_iterations_for_bayesian_fit,
-                         params_initial_fit_bayesian, type_for_MLE, params_class)
+                         params_initial_fit_bayesian, type_for_MLE, params_class, temporal_covariate_for_fit)
         self.max_degree = max_degree
 
     @property
diff --git a/extreme_fit/model/margin_model/polynomial_margin_model/spatio_temporal_polynomial_model.py b/extreme_fit/model/margin_model/polynomial_margin_model/spatio_temporal_polynomial_model.py
index b474c213fc35b2081e702be4fd6a540f2201a1c8..85a6a4000e361a0d788a21da7b1c56b4e8dbd352 100644
--- a/extreme_fit/model/margin_model/polynomial_margin_model/spatio_temporal_polynomial_model.py
+++ b/extreme_fit/model/margin_model/polynomial_margin_model/spatio_temporal_polynomial_model.py
@@ -8,9 +8,10 @@ class AbstractSpatioTemporalPolynomialModel(PolynomialMarginModel):
 
     def __init__(self, coordinates: AbstractCoordinates, params_user=None, starting_point=None,
                  fit_method=MarginFitMethod.extremes_fevd_mle, nb_iterations_for_bayesian_fit=5000,
-                 params_initial_fit_bayesian=None, type_for_MLE="GEV", params_class=GevParams, max_degree=4):
+                 params_initial_fit_bayesian=None, type_for_MLE="GEV", params_class=GevParams, max_degree=4,
+                 temporal_covariate_for_fit=None):
         super().__init__(coordinates, params_user, starting_point, fit_method, nb_iterations_for_bayesian_fit,
-                         params_initial_fit_bayesian, type_for_MLE, params_class, max_degree)
+                         params_initial_fit_bayesian, type_for_MLE, params_class, max_degree, temporal_covariate_for_fit)
         self.drop_duplicates = False
 
 
diff --git a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
index 5508f06c0f8a5fc6a0ef68a20a1ce1e418522835..f34480115283fb87afca6cf768cd46b4d3cfe65c 100644
--- a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
+++ b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py
@@ -6,6 +6,8 @@ from extreme_fit.model.margin_model.polynomial_margin_model.utils import ALTITUD
     ALTITUDINAL_GEV_MODELS_LOCATION_QUADRATIC_MINIMUM, ALTITUDINAL_GEV_MODELS_LOCATION_ONLY_SCALE_ALTITUDES, ALTITUDINAL_GEV_MODELS_LOCATION
 from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.plot_total_aic import plot_total_aic, \
     plot_individual_aic
+from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import \
+    MeanGlobalTemperatureCovariate
 
 mpl.rcParams['text.usetex'] = True
 mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
@@ -42,7 +44,8 @@ def plot_altitudinal_fit(studies, massif_names=None):
     visualizer = AltitudesStudiesVisualizerForNonStationaryModels(studies=studies,
                                                                   model_classes=model_classes,
                                                                   massif_names=massif_names,
-                                                                  show=False)
+                                                                  show=False,
+                                                                  temporal_covariate_for_fit=None)
     # Plot the results for the model that minimizes the individual aic
     plot_individual_aic(visualizer)
     # Plot the results for the model that minimizes the total aic
@@ -55,10 +58,10 @@ def main():
     #  vient probablement de certains zéros pas pris en compte pour le fit,
     # mais pris en compte pour le calcul de mon aic
     # altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:]
-    altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900, 4200, 4500, 4800][:]
     # altitudes = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900][:]
     # altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][:]
     # altitudes = [600, 900, 1200, 1500, 1800, 2100][:]
+    altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900, 4200, 4500, 4800][:]
     # altitudes = [1500, 1800][:]
     study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:2]
     study_classes = [SafranPrecipitation1Day, SafranPrecipitation3Days, SafranPrecipitation5Days,
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
index db0e3634ea867332ec751352701ceaaa75f2bddc..38afd43d980db6ed995b384676cac7962146e3dc 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
@@ -26,11 +26,13 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
                  show=False,
                  massif_names=None,
                  fit_method=MarginFitMethod.extremes_fevd_mle,
+                 temporal_covariate_for_fit=None,
                  display_only_model_that_pass_anderson_test=True):
         super().__init__(studies.study, show=show, save_to_file=not show)
         self.studies = studies
         self.non_stationary_models = model_classes
         self.fit_method = fit_method
+        self.temporal_covariate_for_fit = temporal_covariate_for_fit
         self.display_only_model_that_pass_anderson_test = display_only_model_that_pass_anderson_test
         self.massif_names = massif_names if massif_names is not None else self.study.all_massif_names()
         self.massif_name_to_massif_id = {m: i for i, m in enumerate(self.massif_names)}
@@ -38,7 +40,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
         self._massif_name_to_one_fold_fit = {}
         for massif_name in self.massif_names:
             dataset = studies.spatio_temporal_dataset(massif_name=massif_name)
-            old_fold_fit = OneFoldFit(massif_name, dataset, model_classes, self.fit_method)
+            old_fold_fit = OneFoldFit(massif_name, dataset, model_classes, self.fit_method, self.temporal_covariate_for_fit)
             self._massif_name_to_one_fold_fit[massif_name] = old_fold_fit
 
     @property
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
index 2687b7ce0bf27679bf6fb9d6fca252b2ef35bbe2..3e8b04095d83ebd0ef6a69882f28ccc64b2fbf54 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
@@ -21,18 +21,21 @@ class OneFoldFit(object):
     SIGNIFICANCE_LEVEL = 0.05
     best_estimator_minimizes_total_aic = False
 
-    def __init__(self, massif_name: str, dataset: AbstractDataset, models_classes, fit_method=MarginFitMethod.extremes_fevd_mle):
+    def __init__(self, massif_name: str, dataset: AbstractDataset, models_classes,
+                 fit_method=MarginFitMethod.extremes_fevd_mle, temporal_covariate_for_fit=None):
         self.massif_name = massif_name
         self.dataset = dataset
         self.models_classes = models_classes
         self.fit_method = fit_method
+        self.temporal_covariate_for_fit = temporal_covariate_for_fit
 
         # Fit Estimators
         self.model_class_to_estimator = {}
         for model_class in models_classes:
             self.model_class_to_estimator[model_class] = fitted_linear_margin_estimator_short(model_class=model_class,
                                                                                               dataset=self.dataset,
-                                                                                              fit_method=self.fit_method)
+                                                                                              fit_method=self.fit_method,
+                                                                                              temporal_covariate_for_fit=self.temporal_covariate_for_fit)
 
         # Best estimator definition
         self.best_estimator_class_for_total_aic = None