Repository: sky-uk/anticipy

Issue:
Improve calendar models
Calendar models need further testing and improvements. Some points that need addressing:
- Current predefined calendar models have a separate parameter per holiday. A simpler model would be better. I would explore using one parameter for Christmas, and another for the remaining holidays.
- There may be issues fitting the model with multiplicative composition
- Add Test fitting a calendar model

Discussion:
codecov-io: # [Codecov](https://codecov.io/gh/sky-uk/anticipy/pull/129?src=pr&el=h1) Report
> Merging [#129](https://codecov.io/gh/sky-uk/anticipy/pull/129?src=pr&el=desc) into [master](https://codecov.io/gh/sky-uk/anticipy/commit/cd07d56e9fd5382bf5bb84b21808e1ed2d48ceb4?src=pr&el=desc) will **decrease** coverage by `0.19%`.
> The diff coverage is `100%`.

[![Impacted file tree graph](https://codecov.io/gh/sky-uk/anticipy/pull/129/graphs/tree.svg?width=650&token=QEpdZ0YcB4&height=150&src=pr)](https://codecov.io/gh/sky-uk/anticipy/pull/129?src=pr&el=tree)

```diff
@@            Coverage Diff            @@
##           master     #129     +/-   ##
=========================================
- Coverage   85.41%   85.21%   -0.2%     
=========================================
  Files           6        6             
  Lines        1419     1434     +15     
=========================================
+ Hits         1212     1222     +10     
- Misses        207      212      +5
```


| [Impacted Files](https://codecov.io/gh/sky-uk/anticipy/pull/129?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [anticipy/forecast\_models.py](https://codecov.io/gh/sky-uk/anticipy/pull/129/diff?src=pr&el=tree#diff-YW50aWNpcHkvZm9yZWNhc3RfbW9kZWxzLnB5) | `87.76% <100%> (-0.49%)` | :arrow_down: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/sky-uk/anticipy/pull/129?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/sky-uk/anticipy/pull/129?src=pr&el=footer). Last update [cd07d56...5f1acdd](https://codecov.io/gh/sky-uk/anticipy/pull/129?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).


Patch:
diff --git a/anticipy/forecast.py b/anticipy/forecast.py
index bbd362e..d7e7d5b 100644
--- a/anticipy/forecast.py
+++ b/anticipy/forecast.py
@@ -789,7 +789,8 @@ def run_forecast(df_y, l_model_trend=None, l_model_season=None,
                  l_season_yearly=None,
                  l_season_weekly=None,
                  verbose=None,
-                 l_model_naive=None
+                 l_model_naive=None,
+                 l_model_calendar=None
                  ):
     """
     Generate forecast for one or more input time series
@@ -888,27 +889,30 @@ def run_forecast(df_y, l_model_trend=None, l_model_season=None,
                                    find_outliers,
                                    l_season_yearly,
                                    l_season_weekly,
-                                   l_model_naive=l_model_naive
+                                   l_model_naive=l_model_naive,
+                                   l_model_calendar=l_model_calendar
                                    )
     else:
         for src_tmp in df_y.source.drop_duplicates():
             if verbose:
                 logger.info('Running forecast for source: %s', src_tmp)
             df_y_tmp = df_y.loc[df_y.source == src_tmp].reset_index(drop=True)
-            dict_result_tmp = run_forecast_single(df_y_tmp,
-                                                  l_model_trend,
-                                                  l_model_season,
-                                                  date_start_actuals,
-                                                  src_tmp,
-                                                  extrapolate_years,
-                                                  season_add_mult,
-                                                  include_all_fits,
-                                                  False,  # Simplify output
-                                                  find_outliers,
-                                                  l_season_yearly,
-                                                  l_season_weekly,
-                                                  l_model_naive=l_model_naive
-                                                  )
+            dict_result_tmp = run_forecast_single(
+                df_y_tmp,
+                l_model_trend,
+                l_model_season,
+                date_start_actuals,
+                src_tmp,
+                extrapolate_years,
+                season_add_mult,
+                include_all_fits,
+                False,  # Simplify output
+                find_outliers,
+                l_season_yearly,
+                l_season_weekly,
+                l_model_naive=l_model_naive,
+                l_model_calendar=l_model_calendar
+            )
             l_dict_result += [dict_result_tmp]
     # Generate output
     dict_result = aggregate_forecast_dict_results(l_dict_result)
@@ -969,7 +973,8 @@ def run_forecast_single(df_y,
                         find_outliers=False,
                         l_season_yearly=None,
                         l_season_weekly=None,
-                        l_model_naive=None
+                        l_model_naive=None,
+                        l_model_calendar=None,
                         ):
     """
     Generate forecast for one input time series
@@ -1086,11 +1091,34 @@ def run_forecast_single(df_y,
     l_df_data += [df_actuals]
 
     if l_model_trend is None:
+        # By default use linear trend, and piecewise linear for series > 2y
+        if 'date' in df_y.columns:
+            # Add calendar models
+            s_date_tmp = df_y.date
+            if 'weight' in df_y.columns:
+                s_date_tmp = s_date_tmp.loc[df_y.weight > 0]
+
+            s_date = s_date_tmp.sort_values().drop_duplicates()
+            max_date_delta = s_date.max() - s_date.min()
+
+            if pd.isna(max_date_delta):
+                use_ramp = False
+            else:
+                use_ramp = (
+                    # Need more than 2 full years
+                    (max_date_delta > pd.Timedelta(2 * 365, unit='d'))
+                )
+        else:
+            use_ramp = False
         # By default, try linear and piecewise linear
-        l_model_trend = [
-            # forecast_models.model_naive,
-            forecast_models.model_linear,
-            forecast_models.model_linear + forecast_models.model_ramp]
+        if use_ramp:
+            l_model_trend = [
+                # forecast_models.model_naive,
+                forecast_models.model_linear,
+                forecast_models.model_linear + forecast_models.model_ramp]
+        else:
+            l_model_trend = [forecast_models.model_linear]
+
     l_model_season_add = None
     l_model_season_mult = None
     if l_model_season is None:
@@ -1118,6 +1146,33 @@ def run_forecast_single(df_y,
     l_model_add = get_list_model(l_model_trend, l_model_season_add, 'add')
     l_model_mult = get_list_model(l_model_trend, l_model_season_mult, 'mult')
 
+    if l_model_calendar is not None and 'date' in df_y.columns:
+        # Add calendar models
+
+        s_date_tmp = df_y.date
+        if 'weight' in df_y.columns:
+            s_date_tmp = s_date_tmp.loc[df_y.weight > 0]
+
+        # Requires daily samples, at least one year worth of data
+        # TODO: Should also check for missing calendar dates
+        s_date = s_date_tmp.sort_values().drop_duplicates()
+        min_date_delta = s_date.diff().min()
+        max_date_delta = s_date.max() - s_date.min()
+
+        if pd.isna(min_date_delta) or pd.isna(max_date_delta):
+            use_calendar = False
+        else:
+            use_calendar = (
+                # Need more than a full year
+                (max_date_delta > pd.Timedelta(365, unit='d')) &
+                # Need at least daily samples
+                (min_date_delta <= pd.Timedelta(1, unit='d'))
+            )
+        if use_calendar:
+            l_model_add = get_list_model(l_model_add, l_model_calendar, 'add')
+            l_model_mult = get_list_model(
+                l_model_mult, l_model_calendar, 'mult')
+
     if season_add_mult == 'add':
         l_model = l_model_add
     elif season_add_mult == 'mult':
diff --git a/anticipy/forecast_models.py b/anticipy/forecast_models.py
index fe7313e..6a50fa2 100644
--- a/anticipy/forecast_models.py
+++ b/anticipy/forecast_models.py
@@ -1437,6 +1437,8 @@ def get_f_model_dummy(dummy):
         f_dummy = dummy
     elif isinstance(dummy, Holiday):
         f_dummy = get_f_dummy_from_holiday(dummy)
+    elif isinstance(dummy, AbstractHolidayCalendar):
+        f_dummy = get_f_dummy_from_calendar(dummy)
     else:
         # If dummy is a list, convert to function
         f_dummy = get_f_dummy_from_list(dummy)
@@ -1495,6 +1497,21 @@ def get_f_dummy_from_list(list_check):
                 'date-like values: %s', list_check)
 
 
+def get_f_dummy_from_calendar(calendar):
+    # Generate dummy model function from a pandas HolidayCalendar
+
+    def f_dummy_calendar(a_x, a_date, **kwargs):
+        # TODO: If we can pass dict_cal as an argument,
+        #       use pre-loaded list of dates for performance
+
+        # TODO: If we can guarantee sorted dates,
+        #       change this to a_date[0], a_date[-1] for performance
+        list_check_date = calendar.holidays(a_date.min(), a_date.max())
+        return np.isin(a_date, list_check_date).astype(float)
+
+    return f_dummy_calendar
+
+
 def get_f_dummy_from_holiday(holiday):
     def f_dummy_holiday(a_x, a_date, **kwargs):
         # TODO: If we can pass dict_cal as an argument,
@@ -1526,11 +1543,8 @@ model_dummy_month_start = get_model_dummy(
     (a_date.day == 1).astype(float))
 
 
-class UKCalendar(AbstractHolidayCalendar):
+class CalendarBankHolUK(AbstractHolidayCalendar):
     rules = [
-        Holiday('New Year\'s Day', month=1, day=1, observance=nearest_workday),
-        Holiday('Christmas', month=12, day=25, observance=nearest_workday),
-        Holiday('Boxing Day', month=10, day=26, observance=nearest_workday),
         GoodFriday,
         EasterMonday,
         # Early May Bank Holiday - first Monday in May
@@ -1542,15 +1556,23 @@ class UKCalendar(AbstractHolidayCalendar):
                 offset=DateOffset(weekday=MO(-1))
                 ),
         # August Bank holiday - Last Monday in August
-        Holiday('August Bank Holiday', month=8, day=1,
-                offset=DateOffset(weekday=MO(1))
+        Holiday('August Bank Holiday', month=8, day=30,
+                offset=DateOffset(weekday=MO(-1))
                 )
     ]
 
 
-def get_model_from_calendar(calendar):
+class CalendarChristmasUK(AbstractHolidayCalendar):
+    rules = [
+        Holiday('New Year\'s Day', month=1, day=1, observance=nearest_workday),
+        Holiday('Christmas', month=12, day=25, observance=nearest_workday),
+        Holiday('Boxing Day', month=12, day=26, observance=nearest_workday),
+    ]
+
+
+def get_model_from_calendars(l_calendar, name=None):
     """
-    Create a ForecastModel based on a pandas Calendar.
+    Create a ForecastModel based on a list of pandas Calendars.
 
     :param calendar:
     :type calendar: pandas.tseries.AbstractHolidayCalendar
@@ -1575,9 +1597,20 @@ def get_model_from_calendar(calendar):
         model_calendar = get_model_from_calendar(USFederalHolidayCalendar())
 
     """
-    l_model_dummy = [get_model_dummy(rule.name, rule)
-                     for rule in calendar.rules]
-    assert (len(l_model_dummy)), 'Need 1+ rules in calendar'
+
+    if isinstance(l_calendar, AbstractHolidayCalendar):
+        l_calendar = [l_calendar]
+
+    # Filter out calendars without rules
+    l_calendar = [calendar for calendar in l_calendar if calendar.rules]
+
+    assert len(l_calendar), 'Need 1+ valid calendars'
+
+    if name is None:
+        name = l_calendar[0].name
+
+    l_model_dummy = [get_model_dummy(calendar.name, calendar)
+                     for calendar in l_calendar]
     f_model_prod = np.prod(l_model_dummy)
     f_model_sum = np.sum(l_model_dummy)
 
@@ -1593,7 +1626,7 @@ def get_model_from_calendar(calendar):
         return f_all_dummies(a_x, a_date, params, is_mult, **kwargs)
 
     model_calendar = ForecastModel(
-        calendar.name,
+        name,
         len(l_model_dummy),
         _f_model_calendar,
         _f_init_params_calendar
@@ -1601,8 +1634,10 @@ def get_model_from_calendar(calendar):
     return model_calendar
 
 
-model_ukcalendar = get_model_from_calendar(UKCalendar())
-model_uscalendar = get_model_from_calendar(USFederalHolidayCalendar())
+model_calendar_uk = get_model_from_calendars(
+    [CalendarChristmasUK(), CalendarBankHolUK()], 'calendar_uk')
+model_calendar_us = get_model_from_calendars(USFederalHolidayCalendar(),
+                                             'calendar_us')
 
 
 def get_model_from_datelist(name=None, *args):


Test patch:
diff --git a/tests/test_forecast_model.py b/tests/test_forecast_model.py
index 90f3ebd..d6e0202 100644
--- a/tests/test_forecast_model.py
+++ b/tests/test_forecast_model.py
@@ -258,24 +258,24 @@ class TestForecastModel(PandasTest):
                     is_mult),
                 a)
 
-        test_model('ukcalendar', model_ukcalendar,
-                   array_ones_in_indices(8, 0),
+        test_model('ukcalendar', model_calendar_uk,
+                   [1, 1],
                    # First parameter changes value of New Year
                    array_ones_in_indices(10, 0) + np.ones(10),
                    l_is_mult=[True])
-        test_model('ukcalendar', model_ukcalendar,
-                   array_ones_in_indices(8, 0),
+        test_model('ukcalendar', model_calendar_uk,
+                   [1, 1],
                    # First parameter changes value of New Year
                    array_ones_in_indices(10, 0) + np.zeros(10),
                    l_is_mult=[False])
 
-        test_model('uscalendar', model_uscalendar,
-                   array_ones_in_indices(10, 0),
+        test_model('uscalendar', model_calendar_us,
+                   [1],
                    # First parameter changes value of New Year
                    array_ones_in_indices(10, 0) + np.ones(10),
                    l_is_mult=[True])
-        test_model('uscalendar', model_uscalendar,
-                   array_ones_in_indices(10, 0),
+        test_model('uscalendar', model_calendar_us,
+                   [1],
                    # First parameter changes value of New Year
                    array_ones_in_indices(10, 0) + np.zeros(10),
                    l_is_mult=[False])
@@ -996,9 +996,18 @@ class TestForecastModel(PandasTest):
             model_season_wday.validate_input(
                 None, None, a_date_complete))
 
-    def test_get_model_from_calendar(self):
-        model_calendar = get_model_from_calendar(UKCalendar())
+    def test_get_model_from_calendars(self):
+        model_calendar = get_model_from_calendars(CalendarChristmasUK())
         logger_info('model_calendar:', model_calendar)
+        self.assertEquals(model_calendar.n_params, 1)
+        logger_info('parameters:', model_calendar.n_params)
+
+        model_calendar = get_model_from_calendars(
+            [CalendarChristmasUK(), CalendarBankHolUK()],
+            'calendar2'
+        )
+        logger_info('model_calendar:', model_calendar)
+        self.assertEquals(model_calendar.n_params, 2)
         logger_info('parameters:', model_calendar.n_params)
 
     def test_get_model_from_date_list(self):