diff --git a/mlair/helpers/statistics.py b/mlair/helpers/statistics.py
index a1e713a8c135800d02ff7c27894485a5da7fae37..b63cbe801abd9f0d325a689a959bc494ad1b4485 100644
--- a/mlair/helpers/statistics.py
+++ b/mlair/helpers/statistics.py
@@ -301,7 +301,7 @@ class SkillScores:
                                                             observation_name=self.observation_name)
                                    for (first, second) in combinations]
         return skill_score
-
+    '''
     def climatological_skill_scores(self, internal_data: Data, forecast_name: str) -> xr.DataArray:
         """
         Calculate climatological skill scores according to Murphy (1988).
@@ -341,6 +341,50 @@ class SkillScores:
                     external_data=external_data).values.flatten())
 
         return skill_score
+    '''
+
+    def climatological_skill_scores(self, internal_data: Data, forecast_name: str) -> xr.DataArray:
+        """
+        Calculate climatological skill scores according to Murphy (1988).
+
+        Calculate all CASES I - IV and terms [ABC][I-IV]. Internal data has to be set by initialisation, external data
+        is part of parameters.
+
+        :param internal_data: internal data
+        :param forecast_name: name of the forecast to use for this calculation (must be available in `data`)
+
+        :return: all CASES as well as all terms
+        """
+        if self.external_data is None:
+            ahead_names = []
+        else:
+            ahead_names = list(self.external_data[self.ahead_dim].data)
+
+        all_terms = ['AI', 'AII', 'AIII', 'AIV', 'BI', 'BII', 'BIV', 'CI', 'CIV', 'CASE I', 'CASE II', 'CASE III',
+                     'CASE IV']
+        skill_score = xr.DataArray(np.full((len(all_terms), len(ahead_names)), np.nan), coords=[all_terms, ahead_names],
+                                   dims=['terms', self.ahead_dim])
+
+        for iahead in ahead_names:
+            data = internal_data.sel({self.ahead_dim: iahead})
+
+            skill_score.loc[["CASE I", "AI", "BI", "CI"], iahead] = np.stack(self._climatological_skill_score(
+                data, mu_type=1, forecast_name=forecast_name, observation_name=self.observation_name).values.flatten())
+
+            skill_score.loc[["CASE II", "AII", "BII"], iahead] = np.stack(self._climatological_skill_score(
+                data, mu_type=2, forecast_name=forecast_name, observation_name=self.observation_name).values.flatten())
+
+            if self.external_data is not None and self.observation_name in self.external_data.coords["type"]:
+                external_data = self.external_data.sel({self.ahead_dim: iahead, "type": [self.observation_name]})
+                skill_score.loc[["CASE III", "AIII"], iahead] = np.stack(self._climatological_skill_score(
+                    data, mu_type=3, forecast_name=forecast_name, observation_name=self.observation_name,
+                    external_data=external_data).values.flatten())
+
+                skill_score.loc[["CASE IV", "AIV", "BIV", "CIV"], iahead] = np.stack(self._climatological_skill_score(
+                    data, mu_type=4, forecast_name=forecast_name, observation_name=self.observation_name,
+                    external_data=external_data).values.flatten())
+
+        return skill_score
 
     def _climatological_skill_score(self, internal_data, observation_name, forecast_name, mu_type=1,
                                     external_data=None):