diff --git a/experiment/paper1_steps/poster_EVAN2019/some_experiment_EVAN.py b/experiment/paper1_steps/poster_EVAN2019/some_experiment_EVAN.py
index e4c731a1952208f6439d598b133472a3279c60cc..6ea0da2462482342c083b7aed7611877034756d5 100644
--- a/experiment/paper1_steps/poster_EVAN2019/some_experiment_EVAN.py
+++ b/experiment/paper1_steps/poster_EVAN2019/some_experiment_EVAN.py
@@ -18,17 +18,19 @@ mpl.rcParams['hatch.linewidth'] = 0.3
 
 
 def main_non_stationary_model_comparison():
-    stop_loop = False
+    stop_loop = True
     for altitude in POSTER_ALTITUDES[:]:
-        for trend_test_class in [GevLocationTrendTest, GevScaleTrendTest, GevLocationAndScaleTrendTest, ComparisonAgainstMu, ComparisonAgainstSigma][::-1]:
+        for trend_test_class in [GevLocationTrendTest, GevScaleTrendTest, GevLocationAndScaleTrendTest,
+                                 ComparisonAgainstMu, ComparisonAgainstSigma][::-1][-1:]:
             vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
                                                       exact_starting_year=1958, reduce_strength_array=False,
                                                       trend_test_class=trend_test_class,
                                                       )
-            # vizualiser.save_to_file = False
+            vizualiser.save_to_file = False
             vizualiser.visualize_massif_trend_test_one_altitude(poster_plot=True, write_text_on_massif=False)
             if stop_loop:
                 return
 
+
 if __name__ == '__main__':
     main_non_stationary_model_comparison()
diff --git a/experiment/trend_analysis/univariate_test/abstract_comparison_non_stationary_model.py b/experiment/trend_analysis/univariate_test/abstract_comparison_non_stationary_model.py
index cb569655c2e2cf5b918d83a1755148ad946a3735..c9fb5671282201d8e48a28623d4db4897fd23b22 100644
--- a/experiment/trend_analysis/univariate_test/abstract_comparison_non_stationary_model.py
+++ b/experiment/trend_analysis/univariate_test/abstract_comparison_non_stationary_model.py
@@ -2,6 +2,7 @@ from experiment.trend_analysis.univariate_test.abstract_gev_trend_test import Ab
 from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters import GevLocationAndScaleTrendTest
 from extreme_estimator.extreme_models.margin_model.temporal_linear_margin_model import \
     NonStationaryLocationAndScaleModel, NonStationaryLocationStationModel, NonStationaryScaleStationModel
+import numpy as np
 
 
 class AbstractComparisonNonStationaryModel(AbstractGevTrendTest):
@@ -17,6 +18,12 @@ class AbstractComparisonNonStationaryModelOneParameter(AbstractComparisonNonStat
     def degree_freedom_chi2(self) -> int:
         return 1
 
+    @property
+    def test_sign(self) -> int:
+        # Test sign correspond to the difference between the 2 likelihoods
+        # Therefore, colors sum up which non stationary model explain best the data
+        return np.sign(self.likelihood_ratio)
+
 
 class ComparisonAgainstMu(AbstractComparisonNonStationaryModelOneParameter, GevLocationAndScaleTrendTest):
 
diff --git a/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py b/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py
index deae22fc07f4e43c8b80bf2f8338279e3ea728d6..95a4345c8ec09ad0fe76c3d11098661632666f99 100644
--- a/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py
+++ b/experiment/trend_analysis/univariate_test/abstract_gev_trend_test.py
@@ -1,5 +1,6 @@
 import numpy as np
 import pandas as pd
+from cached_property import cached_property
 from scipy.stats import chi2
 
 from experiment.trend_analysis.univariate_test.abstract_univariate_test import AbstractUnivariateTest
@@ -109,11 +110,11 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
     def test_sign(self) -> int:
         return np.sign(self.test_trend_slope_strength)
 
-    def get_non_stationary_linear_coef(self, gev_param_name):
+    def get_non_stationary_linear_coef(self, gev_param_name: str):
         return self.non_stationary_estimator.margin_function_fitted.get_coef(gev_param_name,
                                                                              AbstractCoordinates.COORDINATE_T)
 
-    @property
+    @cached_property
     def non_stationary_constant_gev_params(self) -> GevParams:
         # Constant parameters correspond to the gev params in 1958
         return self.non_stationary_estimator.margin_function_fitted.get_gev_params(coordinate=np.array([1958]),
@@ -145,8 +146,8 @@ class AbstractGevTrendTest(AbstractUnivariateTest):
     def variance_difference_same_sign_as_slope_strenght(self) -> bool:
         return False
 
-    def mean_difference(self, zeta0, mu1=0.0, sigma1=0.0):
-        return GevParams(loc=mu1, scale=sigma1, shape=zeta0).mean
+    def mean_difference(self, zeta0: float, mu1: float = 0.0, sigma1: float = 0.0) -> float:
+        return GevParams(loc=mu1, scale=sigma1, shape=zeta0, accept_zero_scale_parameter=True).mean
 
     @property
     def test_trend_constant_quantile(self):
diff --git a/extreme_estimator/margin_fits/extreme_params.py b/extreme_estimator/margin_fits/extreme_params.py
index a20ace4f1493ebe88aa2ab28b29524b8fae33552..026c40405480753082e75b99b848e6a2e0864f21 100644
--- a/extreme_estimator/margin_fits/extreme_params.py
+++ b/extreme_estimator/margin_fits/extreme_params.py
@@ -14,7 +14,7 @@ class ExtremeParams(AbstractParams, ABC):
         self.location = loc
         self.scale = scale
         self.shape = shape
-        # Scale cannot be negative
+        # By default, scale cannot be negative
         # (sometimes it happens, when we want to find a quantile for every point of a 2D map
         # then it can happen that a corner point that was not used for fitting correspond to a negative scale,
         # in the case we set all the parameters as equal to np.nan, and we will not display those points)
diff --git a/extreme_estimator/margin_fits/gev/gev_params.py b/extreme_estimator/margin_fits/gev/gev_params.py
index 769df9fdbef859cf5f3fa552dfa30d01ef335195..395032b551a9f18ff10ae347b2f9d03465ae705c 100644
--- a/extreme_estimator/margin_fits/gev/gev_params.py
+++ b/extreme_estimator/margin_fits/gev/gev_params.py
@@ -1,4 +1,5 @@
 from collections import OrderedDict
+from typing import List
 
 from cached_property import cached_property
 
@@ -15,9 +16,11 @@ class GevParams(ExtremeParams):
     SUMMARY_NAMES = PARAM_NAMES + ExtremeParams.QUANTILE_NAMES
     NB_SUMMARY_NAMES = len(SUMMARY_NAMES)
 
-    def __init__(self, loc: float, scale: float, shape: float, block_size: int = None):
+    def __init__(self, loc: float, scale: float, shape: float, block_size: int = None, accept_zero_scale_parameter=False):
         super().__init__(loc, scale, shape)
         self.block_size = block_size
+        if accept_zero_scale_parameter and scale == 0.0:
+            self.has_undefined_parameters = False
 
     def quantile(self, p) -> float:
         if self.has_undefined_parameters:
@@ -46,6 +49,7 @@ class GevParams(ExtremeParams):
         return self.to_dict().__str__()
 
     def quantile_strength_evolution(self, p=0.99, mu1=0.0, sigma1=0.0):
+        # todo: chagne the name to time derivative
         """
         Compute the variation of some quantile with respect to time.
         (when mu1 and sigma1 can be modified with time)
@@ -63,13 +67,13 @@ class GevParams(ExtremeParams):
 
     # Compute some indicators (such as the mean and the variance)
 
-    def g(self, k):
+    def g(self, k) -> float:
         # Compute the g_k parameters as defined in wikipedia
         # https://fr.wikipedia.org/wiki/Loi_d%27extremum_g%C3%A9n%C3%A9ralis%C3%A9e
         return gamma(1 - k * self.shape)
 
     @property
-    def mean(self):
+    def mean(self) -> float:
         if self.has_undefined_parameters:
             return np.nan
         elif self.shape >= 1:
@@ -78,7 +82,7 @@ class GevParams(ExtremeParams):
             return self.location + self.scale * (self.g(k=1) - 1) / self.shape
 
     @property
-    def variance(self):
+    def variance(self) -> float:
         if self.has_undefined_parameters:
             return np.nan
         elif self.shape >= 0.5:
@@ -87,11 +91,11 @@ class GevParams(ExtremeParams):
             return ((self.scale / self.shape) ** 2) * (self.g(k=2) - self.g(k=1) ** 2)
 
     @property
-    def std(self):
+    def std(self) -> float:
         return np.sqrt(self.variance)
 
     @classmethod
-    def indicator_names(cls):
+    def indicator_names(cls) -> List[str]:
         return ['mean', 'std'] + cls.QUANTILE_NAMES[:2]
 
     @cached_property