bshishov · July 22, 2025 10:05 · Sep 2, 2023 · Nov 12, 2022 · Apr 25, 2018
diff --git a/forecasting_metrics.py b/forecasting_metrics.py
@@ -93,7 +93,10 @@ def mae(actual: np.ndarray, predicted: np.ndarray):
     return np.mean(np.abs(_error(actual, predicted)))
 
 
-mad = mae  # Mean Absolute Deviation (it is the same as MAE)
+def mad(actual: np.ndarray, predicted: np.ndarray):
+    """ Mean Absolute Deviation """
+    error = _error(actual, predicted)
+    return np.mean(np.abs(error - np.mean(error)))
 
 
 def gmae(actual: np.ndarray, predicted: np.ndarray):

diff --git a/forecasting_metrics.py b/forecasting_metrics.py
@@ -255,7 +255,7 @@ def umbrae(actual: np.ndarray, predicted: np.ndarray, benchmark: np.ndarray = No
 
 def mda(actual: np.ndarray, predicted: np.ndarray):
     """ Mean Directional Accuracy """
-    return np.mean((np.sign(actual[1:] - actual[:-1]) == np.sign(predicted[1:] - predicted[:-1])).astype(int))
+    return np.mean((np.sign(actual[1:] - actual[:-1]) == np.sign(predicted[1:] - actual[:-1])).astype(int))
 
 
 METRICS = {

diff --git a/forecasting_metrics.py b/forecasting_metrics.py
@@ -0,0 +1,307 @@
+import numpy as np
+
+EPSILON = 1e-10
+
+
+def _error(actual: np.ndarray, predicted: np.ndarray):
+    """ Simple error """
+    return actual - predicted
+
+
+def _percentage_error(actual: np.ndarray, predicted: np.ndarray):
+    """
+    Percentage error
+
+    Note: result is NOT multiplied by 100
+    """
+    return _error(actual, predicted) / (actual + EPSILON)
+
+
+def _naive_forecasting(actual: np.ndarray, seasonality: int = 1):
+    """ Naive forecasting method which just repeats previous samples """
+    return actual[:-seasonality]
+
+
+def _relative_error(actual: np.ndarray, predicted: np.ndarray, benchmark: np.ndarray = None):
+    """ Relative Error """
+    if benchmark is None or isinstance(benchmark, int):
+        # If no benchmark prediction provided - use naive forecasting
+        if not isinstance(benchmark, int):
+            seasonality = 1
+        else:
+            seasonality = benchmark
+        return _error(actual[seasonality:], predicted[seasonality:]) /\
+               (_error(actual[seasonality:], _naive_forecasting(actual, seasonality)) + EPSILON)
+
+    return _error(actual, predicted) / (_error(actual, benchmark) + EPSILON)
+
+
+def _bounded_relative_error(actual: np.ndarray, predicted: np.ndarray, benchmark: np.ndarray = None):
+    """ Bounded Relative Error """
+    if benchmark is None or isinstance(benchmark, int):
+        # If no benchmark prediction provided - use naive forecasting
+        if not isinstance(benchmark, int):
+            seasonality = 1
+        else:
+            seasonality = benchmark
+
+        abs_err = np.abs(_error(actual[seasonality:], predicted[seasonality:]))
+        abs_err_bench = np.abs(_error(actual[seasonality:], _naive_forecasting(actual, seasonality)))
+    else:
+        abs_err = np.abs(_error(actual, predicted))
+        abs_err_bench = np.abs(_error(actual, benchmark))
+
+    return abs_err / (abs_err + abs_err_bench + EPSILON)
+
+
+def _geometric_mean(a, axis=0, dtype=None):
+    """ Geometric mean """
+    if not isinstance(a, np.ndarray):  # if not an ndarray object attempt to convert it
+        log_a = np.log(np.array(a, dtype=dtype))
+    elif dtype:  # Must change the default dtype allowing array type
+        if isinstance(a, np.ma.MaskedArray):
+            log_a = np.log(np.ma.asarray(a, dtype=dtype))
+        else:
+            log_a = np.log(np.asarray(a, dtype=dtype))
+    else:
+        log_a = np.log(a)
+    return np.exp(log_a.mean(axis=axis))
+
+
+def mse(actual: np.ndarray, predicted: np.ndarray):
+    """ Mean Squared Error """
+    return np.mean(np.square(_error(actual, predicted)))
+
+
+def rmse(actual: np.ndarray, predicted: np.ndarray):
+    """ Root Mean Squared Error """
+    return np.sqrt(mse(actual, predicted))
+
+
+def nrmse(actual: np.ndarray, predicted: np.ndarray):
+    """ Normalized Root Mean Squared Error """
+    return rmse(actual, predicted) / (actual.max() - actual.min())
+
+
+def me(actual: np.ndarray, predicted: np.ndarray):
+    """ Mean Error """
+    return np.mean(_error(actual, predicted))
+
+
+def mae(actual: np.ndarray, predicted: np.ndarray):
+    """ Mean Absolute Error """
+    return np.mean(np.abs(_error(actual, predicted)))
+
+
+mad = mae  # Mean Absolute Deviation (it is the same as MAE)
+
+
+def gmae(actual: np.ndarray, predicted: np.ndarray):
+    """ Geometric Mean Absolute Error """
+    return _geometric_mean(np.abs(_error(actual, predicted)))
+
+
+def mdae(actual: np.ndarray, predicted: np.ndarray):
+    """ Median Absolute Error """
+    return np.median(np.abs(_error(actual, predicted)))
+
+
+def mpe(actual: np.ndarray, predicted: np.ndarray):
+    """ Mean Percentage Error """
+    return np.mean(_percentage_error(actual, predicted))
+
+
+def mape(actual: np.ndarray, predicted: np.ndarray):
+    """
+    Mean Absolute Percentage Error
+
+    Properties:
+        + Easy to interpret
+        + Scale independent
+        - Biased, not symmetric
+        - Undefined when actual[t] == 0
+
+    Note: result is NOT multiplied by 100
+    """
+    return np.mean(np.abs(_percentage_error(actual, predicted)))
+
+
+def mdape(actual: np.ndarray, predicted: np.ndarray):
+    """
+    Median Absolute Percentage Error
+
+    Note: result is NOT multiplied by 100
+    """
+    return np.median(np.abs(_percentage_error(actual, predicted)))
+
+
+def smape(actual: np.ndarray, predicted: np.ndarray):
+    """
+    Symmetric Mean Absolute Percentage Error
+
+    Note: result is NOT multiplied by 100
+    """
+    return np.mean(2.0 * np.abs(actual - predicted) / ((np.abs(actual) + np.abs(predicted)) + EPSILON))
+
+
+def smdape(actual: np.ndarray, predicted: np.ndarray):
+    """
+    Symmetric Median Absolute Percentage Error
+
+    Note: result is NOT multiplied by 100
+    """
+    return np.median(2.0 * np.abs(actual - predicted) / ((np.abs(actual) + np.abs(predicted)) + EPSILON))
+
+
+def maape(actual: np.ndarray, predicted: np.ndarray):
+    """
+    Mean Arctangent Absolute Percentage Error
+
+    Note: result is NOT multiplied by 100
+    """
+    return np.mean(np.arctan(np.abs((actual - predicted) / (actual + EPSILON))))
+
+
+def mase(actual: np.ndarray, predicted: np.ndarray, seasonality: int = 1):
+    """
+    Mean Absolute Scaled Error
+
+    Baseline (benchmark) is computed with naive forecasting (shifted by @seasonality)
+    """
+    return mae(actual, predicted) / mae(actual[seasonality:], _naive_forecasting(actual, seasonality))
+
+
+def std_ae(actual: np.ndarray, predicted: np.ndarray):
+    """ Normalized Absolute Error """
+    __mae = mae(actual, predicted)
+    return np.sqrt(np.sum(np.square(_error(actual, predicted) - __mae))/(len(actual) - 1))
+
+
+def std_ape(actual: np.ndarray, predicted: np.ndarray):
+    """ Normalized Absolute Percentage Error """
+    __mape = mape(actual, predicted)
+    return np.sqrt(np.sum(np.square(_percentage_error(actual, predicted) - __mape))/(len(actual) - 1))
+
+
+def rmspe(actual: np.ndarray, predicted: np.ndarray):
+    """
+    Root Mean Squared Percentage Error
+
+    Note: result is NOT multiplied by 100
+    """
+    return np.sqrt(np.mean(np.square(_percentage_error(actual, predicted))))
+
+
+def rmdspe(actual: np.ndarray, predicted: np.ndarray):
+    """
+    Root Median Squared Percentage Error
+
+    Note: result is NOT multiplied by 100
+    """
+    return np.sqrt(np.median(np.square(_percentage_error(actual, predicted))))
+
+
+def rmsse(actual: np.ndarray, predicted: np.ndarray, seasonality: int = 1):
+    """ Root Mean Squared Scaled Error """
+    q = np.abs(_error(actual, predicted)) / mae(actual[seasonality:], _naive_forecasting(actual, seasonality))
+    return np.sqrt(np.mean(np.square(q)))
+
+
+def inrse(actual: np.ndarray, predicted: np.ndarray):
+    """ Integral Normalized Root Squared Error """
+    return np.sqrt(np.sum(np.square(_error(actual, predicted))) / np.sum(np.square(actual - np.mean(actual))))
+
+
+def rrse(actual: np.ndarray, predicted: np.ndarray):
+    """ Root Relative Squared Error """
+    return np.sqrt(np.sum(np.square(actual - predicted)) / np.sum(np.square(actual - np.mean(actual))))
+
+
+def mre(actual: np.ndarray, predicted: np.ndarray, benchmark: np.ndarray = None):
+    """ Mean Relative Error """
+    return np.mean(_relative_error(actual, predicted, benchmark))
+
+
+def rae(actual: np.ndarray, predicted: np.ndarray):
+    """ Relative Absolute Error (aka Approximation Error) """
+    return np.sum(np.abs(actual - predicted)) / (np.sum(np.abs(actual - np.mean(actual))) + EPSILON)
+
+
+def mrae(actual: np.ndarray, predicted: np.ndarray, benchmark: np.ndarray = None):
+    """ Mean Relative Absolute Error """
+    return np.mean(np.abs(_relative_error(actual, predicted, benchmark)))
+
+
+def mdrae(actual: np.ndarray, predicted: np.ndarray, benchmark: np.ndarray = None):
+    """ Median Relative Absolute Error """
+    return np.median(np.abs(_relative_error(actual, predicted, benchmark)))
+
+
+def gmrae(actual: np.ndarray, predicted: np.ndarray, benchmark: np.ndarray = None):
+    """ Geometric Mean Relative Absolute Error """
+    return _geometric_mean(np.abs(_relative_error(actual, predicted, benchmark)))
+
+
+def mbrae(actual: np.ndarray, predicted: np.ndarray, benchmark: np.ndarray = None):
+    """ Mean Bounded Relative Absolute Error """
+    return np.mean(_bounded_relative_error(actual, predicted, benchmark))
+
+
+def umbrae(actual: np.ndarray, predicted: np.ndarray, benchmark: np.ndarray = None):
+    """ Unscaled Mean Bounded Relative Absolute Error """
+    __mbrae = mbrae(actual, predicted, benchmark)
+    return __mbrae / (1 - __mbrae)
+
+
+def mda(actual: np.ndarray, predicted: np.ndarray):
+    """ Mean Directional Accuracy """
+    return np.mean((np.sign(actual[1:] - actual[:-1]) == np.sign(predicted[1:] - predicted[:-1])).astype(int))
+
+
+METRICS = {
+    'mse': mse,
+    'rmse': rmse,
+    'nrmse': nrmse,
+    'me': me,
+    'mae': mae,
+    'mad': mad,
+    'gmae': gmae,
+    'mdae': mdae,
+    'mpe': mpe,
+    'mape': mape,
+    'mdape': mdape,
+    'smape': smape,
+    'smdape': smdape,
+    'maape': maape,
+    'mase': mase,
+    'std_ae': std_ae,
+    'std_ape': std_ape,
+    'rmspe': rmspe,
+    'rmdspe': rmdspe,
+    'rmsse': rmsse,
+    'inrse': inrse,
+    'rrse': rrse,
+    'mre': mre,
+    'rae': rae,
+    'mrae': mrae,
+    'mdrae': mdrae,
+    'gmrae': gmrae,
+    'mbrae': mbrae,
+    'umbrae': umbrae,
+    'mda': mda,
+}
+
+
+def evaluate(actual: np.ndarray, predicted: np.ndarray, metrics=('mae', 'mse', 'smape', 'umbrae')):
+    results = {}
+    for name in metrics:
+        try:
+            results[name] = METRICS[name](actual, predicted)
+        except Exception as err:
+            results[name] = np.nan
+            print('Unable to compute metric {0}: {1}'.format(name, err))
+    return results
+
+
+def evaluate_all(actual: np.ndarray, predicted: np.ndarray):
+    return evaluate(actual, predicted, metrics=set(METRICS.keys()))