Repository: TomMonks/forecast-tools (language: python)

Issue:
auto_naive for v0.2
Added auto_naive function to model_selection
Unit tests in tests/test_model_selection
Examples in examples/04_auto_naive.ipynb
Increase unit test coverage > 90%.

Updated readme.

Big fixes:
* MASE when scaled by Naive.



Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
index 0000000..6dd00ea
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,18 @@
+# Change log
+
+## v0.1.5
+
+auto_naive
+increased unit-test coverage (> 90%)
+
+bug fix: mase for seasonal naive
+
+## v0.1.4
+
+cross-validation
+mase
+
+
+
+
+
diff --git a/README.md b/README.md
index 0780cd4..0e27ddf 100644
--- a/README.md
+++ b/README.md
@@ -1,11 +1,51 @@
-# forecast_tools
-
+# forecast-tools: fundamental tools to support the forecasting process in python.
 
 [![DOI](https://zenodo.org/badge/250494795.svg)](https://zenodo.org/badge/latestdoi/250494795)
-
+[![PyPI version fury.io](https://badge.fury.io/py/forecast-tools.svg)](https://pypi.python.org/pypi/forecast-tools/)
+[![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT)
 [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/TomMonks/forecast-tools/master)
+[![Python 3.6+](https://img.shields.io/badge/python-3.6+-blue.svg)](https://www.python.org/downloads/release/python-360+/)
+
+ forecast-tools has been developed to support forecasting education and applied forecasting research.  It is MIT licensed and freely available to practitioners, students and researchers via PyPi.  There is a long term plan to make forecast-tools available via conda-forge.
+
+ # Vision for forecast-tools
+
+ 1. Deliverhigh quality reliable code for forecasting education and practice with full documentation and unit testing.
+ 2. Provide a simple to use pythonic interface that users of `statsmodels` and `sklearn` will recognise.
+ 3. To improve the quality of Machine Learning time series forecasting and encourage the use of best practice.
+
+# Features:
+
+1. Implementation of classic naive forecast benchmarks such as Naive Forecast 1 along with prediction intervals
+2. Implementation of scale-dependent, relative and scaled forecast errors.
+3. Rolling forecast origin and sliding window for time series cross validation
+4. Built in daily level datasets
+
+## Two simple ways to explore forecast-tools
 
-forecast_tools provides fundamental tools to support the forecasting process in python
+1. `pip install forecast-tools`
+2. Click on the launch-binder at the top of this readme. This will open example Jupyter notebooks in the cloud via Binder.
+
+## Citation
+
+If you use forecast-tools for research, a practical report, education or any reason please include the following citation.
+
+> Monks, Thomas. (2020). forecast-tools: fundamental tools to support the forecasting process in python. Zenodo. http://doi.org/10.5281/zenodo.3969789
+
+```tex
+@software{forecast_tools_3969789,
+  author       = {Thomas Monks},
+  title        = {forecast-tools: fundamental tools to support the forecasting process in python},
+  year         = 2020,
+  publisher    = {Zenodo},
+  doi          = {10.5281/zenodo.3969789},
+  url          = {https://doi.org/10.5281/zenodo.3969789}
+}
+```
+
+## Contributing to forecast-tools
+
+Please fork Dev, make your modifications, run the unit tests and submit a pull request for review.
 
 Development environment:
 
@@ -13,6 +53,8 @@ Development environment:
 
 * `conda activate forecast_dev`
 
-Launch in the cloud via Binder.
+Unit tests are provided and can be run from the command `pytest` and its coverage extension.  Run the following in the terminal.
+
+* `pytest --cov=forecast_tools tests/`
 
-* Click on the binder badge at the top of this readme to run the example Jupyter Notebooks in the cloud.
+**All contributions are welcome and must include unit tests!**
\ No newline at end of file
diff --git a/binder/environment.yml b/binder/environment.yml
index c167f01..812ab52 100644
--- a/binder/environment.yml
+++ b/binder/environment.yml
@@ -13,10 +13,10 @@ dependencies:
   - scipy=1.4.1
   - seaborn=0.10.0
   - spyder=4.0.1
-  - statsmodels=0.11.0
+  - statsmodels=0.11.1
   - pip:
     - sklearn==0.0
     - joblib==0.14.1
     - pytest-cov==2.10.0
-    - forecast-tools==0.1.3
+    - forecast-tools==0.1.4
 
diff --git a/examples/01_example_benchmark_forecasts.ipynb b/examples/01_benchmark_forecasts.ipynb
similarity index 100%
rename from examples/01_example_benchmark_forecasts.ipynb
rename to examples/01_benchmark_forecasts.ipynb
diff --git a/examples/02_example_forecast_error.ipynb b/examples/02_forecast_error.ipynb
similarity index 100%
rename from examples/02_example_forecast_error.ipynb
rename to examples/02_forecast_error.ipynb
diff --git a/examples/03_example_cross_validation.ipynb b/examples/03_cross_validation.ipynb
similarity index 100%
rename from examples/03_example_cross_validation.ipynb
rename to examples/03_cross_validation.ipynb
diff --git a/examples/04_auto_naive.ipynb b/examples/04_auto_naive.ipynb
new file mode 100644
index 0000000..29922a9
--- /dev/null
+++ b/examples/04_auto_naive.ipynb
@@ -0,0 +1,301 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Automatically selecting a naive model to use as a benchmark\n",
+    "\n",
+    "forecast-tools provides a `auto_naive` function that uses point-forecast cross validation to select the 'best' naive model to use as a benchmark.  \n",
+    "\n",
+    "The function tests all of the naive `Forecast` methods."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Imports"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from forecast_tools.datasets import load_emergency_dept\n",
+    "from forecast_tools.model_selection import auto_naive                                    "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Help on function auto_naive in module forecast_tools.model_selection:\n",
+      "\n",
+      "auto_naive(y_train, horizon=1, seasonal_period=1, min_train_size='auto', method='cv', step=1, window_size='auto', metric='mae')\n",
+      "    Automatic selection of the 'best' naive benchmark on a 'single' series\n",
+      "    \n",
+      "    The selection process uses out-of-sample cv performance.\n",
+      "    \n",
+      "    By default auto_naive uses cross validation to estimate the mean\n",
+      "    point forecast peformance of all naive methods.  It selects the method\n",
+      "    with the lowest point forecast metric on average.\n",
+      "    \n",
+      "    If there is limited data for training a basic holdout sample could be\n",
+      "    used.\n",
+      "    \n",
+      "    Dev note: the plan is to update this to work with multiple series.\n",
+      "    It would be best to use MASE for multiple series comparison.\n",
+      "    \n",
+      "    Parameters:\n",
+      "    ----------\n",
+      "    y_train: array-like\n",
+      "        training data.  typically in a pandas.Series, pandas.DataFrame\n",
+      "        or numpy.ndarray format.\n",
+      "    \n",
+      "    horizon: int, optional (default=1)\n",
+      "        Forecast horizon.\n",
+      "    \n",
+      "    seasonal_period: int, optional (default=1)\n",
+      "        Frequency of the data.  E.g. 7 for weekly pattern, 12 for monthly\n",
+      "        365 for daily.\n",
+      "    \n",
+      "    min_train_size: int or str, optional (default='auto')\n",
+      "        The size of the initial training set (if method=='ro' or 'sw').\n",
+      "        If 'auto' then then min_train_size is set to len(y_train) // 3\n",
+      "        If main_train_size='auto' and method='holdout' then\n",
+      "        min_train_size = len(y_train) - horizon.\n",
+      "    \n",
+      "    method: str, optional (default='cv')\n",
+      "        out of sample selection method.\n",
+      "        'ro' - rolling forecast origin\n",
+      "        'sw' - sliding window\n",
+      "        'cv' - scores from both ro and sw\n",
+      "        'holdout' - single train/test split\n",
+      "         Methods'ro' and 'sw' are similar, however, sw has a fixed\n",
+      "         window_size and drops older data from training.\n",
+      "    \n",
+      "    step: int, optional (default=1)\n",
+      "        The stride/step of the cross-validation. I.e. the number\n",
+      "        of observations to move forward between folds.\n",
+      "    \n",
+      "    window_size: str or int, optional (default='auto')\n",
+      "        The window_size if using sliding window cross validation\n",
+      "        When 'auto' and method='sw' then\n",
+      "        window_size=len(y_train) // 3\n",
+      "    \n",
+      "    metric: str, optional (default='mae')\n",
+      "        The metric to measure out of sample accuracy.\n",
+      "        Options: mase, mae, mape, smape, mse, rmse, me.\n",
+      "    \n",
+      "    Returns:\n",
+      "    --------\n",
+      "    dict\n",
+      "        'model': baseline.Forecast\n",
+      "        f'{metric}': float\n",
+      "    \n",
+      "        Contains the model and its CV performance.\n",
+      "    \n",
+      "    Raises:\n",
+      "    -------\n",
+      "    ValueError\n",
+      "        For invalid method, metric, window_size parameters\n",
+      "    \n",
+      "    See Also:\n",
+      "    --------\n",
+      "    forecast_tools.baseline.Naive1\n",
+      "    forecast_tools.baseline.SNaive\n",
+      "    forecast_tools.baseline.Drift\n",
+      "    forecast_tools.baseline.Average\n",
+      "    forecast_tools.baseline.EnsembleNaive\n",
+      "    forecast_tools.baseline.baseline_estimators\n",
+      "    forecast_tools.model_selection.rolling_forecast_origin\n",
+      "    forecast_tools.model_selection.sliding_window\n",
+      "    forecast_tools.model_selection.mase_cross_validation_score\n",
+      "    forecast_tools.metrics.mean_absolute_scaled_error\n",
+      "    \n",
+      "    Examples:\n",
+      "    ---------\n",
+      "    Measuring MAE and taking the best method using both\n",
+      "    rolling origin and sliding window cross validation\n",
+      "    of a 56 day forecast.\n",
+      "    \n",
+      "    >>> from forecast_tools.datasets import load_emergency_dept\n",
+      "    >>> y_train = load_emergency_dept\n",
+      "    >>> best = auto_naive(y_train, seasonal_period=7, horizon=56)\n",
+      "    >>> best\n",
+      "    {'model': Average(), 'mae': 19.63791579700355}\n",
+      "    \n",
+      "    \n",
+      "    Take a step of 7 days between cv folds.\n",
+      "    \n",
+      "    >>> from forecast_tools.datasets import load_emergency_dept\n",
+      "    >>> y_train = load_emergency_dept\n",
+      "    >>> best = auto_naive(y_train, seasonal_period=7, horizon=56,\n",
+      "        ...               step=7)\n",
+      "    >>> best\n",
+      "    {'model': Average(), 'mae': 19.675635558539383}\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "help(auto_naive)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Load the training data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "y_train = load_emergency_dept()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Select the best naive model for a h-step horizon of 7 days.\n",
+    "\n",
+    "Let's select a method for the emergency deparment daily level to predict 7 days ahead.  By default the function using the **mean absolute error** to evaluate forecast accuracy."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{'model': Average(), 'mae': 19.679856211931035}"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "best = auto_naive(y_train, horizon=7, seasonal_period=7)\n",
+    "best"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([221.06395349, 221.06395349, 221.06395349, 221.06395349,\n",
+       "       221.06395349, 221.06395349, 221.06395349])"
+      ]
+     },
+     "execution_count": 5,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "y_preds = best['model'].fit_predict(y_train, horizon=7)\n",
+    "y_preds"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Using a different forecasting error metric"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{'model': Average(), 'mape': 8.69955926909263}"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "best = auto_naive(y_train, horizon=7, seasonal_period=7, metric='mape')\n",
+    "best"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Using a single train-test split when data are limited.\n",
+    "\n",
+    "If your forecast horizon means that h-step cross-validation is infeasible then you can automatically select using a single holdout sample."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{'model': Average(), 'mae': 30.182280627384486}"
+      ]
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "best = auto_naive(y_train, horizon=7, seasonal_period=7, method='holdout')\n",
+    "best"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.1"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/forecast_tools/baseline.py b/forecast_tools/baseline.py
index 9a1a08d..6f3f44a 100644
--- a/forecast_tools/baseline.py
+++ b/forecast_tools/baseline.py
@@ -1,5 +1,5 @@
 '''
-Contains classes for creating simple baseline/benchmark forecasts
+Tools for simple baseline/benchmark forecasts
 
 These methods might serve as the forecast themselves, but are more likely
 to be used as a baseline to evaluate if more complex models offer a sufficient
@@ -16,6 +16,9 @@
 
 Drift:
     Carry forward last time period, but allow for upwards/downwards drift.
+
+EnsembleNaive:
+    An unweighted average of all of the Naive forecasting methods.
 '''
 
 import numpy as np
@@ -23,7 +26,6 @@
 from scipy.stats import norm, t
 from abc import ABC, abstractmethod
 
-
 # Boolean, unsigned integer, signed integer, float, complex.
 _NUMERIC_KINDS = set('buifc')
 
@@ -219,15 +221,89 @@ def _std_h(self, horizon):
 
 class Naive1(Forecast):
     '''
-    Naive1 or NF1: Carry the last value foreward across a
+    Naive forecast 1 or NF1: Carry the last value foreward across a
     forecast horizon
 
-    See Makridakis, Wheelwright and Hyndman (1998) and
-    Hyndman:
-    https://otexts.com/fpp2/simple-methods.html
+    For details and theory see [1]
 
-    '''
+    Attributes
+    ----------
+    fittedvalues: pd.Series
+        In-sample predictions of training data
+    resid: pd.Series
+        In-sample residuals
+
+    Methods
+    -------
+    fit(train)
+        fit the model to training data
+    predict(horizon, return_predict_int=False, alpha=None)
+        Predict h-steps ahead
+    fit_predict(train, horizons, return_predict_int=False, alpha=None)
+        convenience method.  combine fit() and predict()
+
+    See Also
+    --------
+    forecast_tools.baseline.SNaive
+    forecast_tools.baseline.Drift
+    forecast_tools.baseline.Average
+    forecast_tools.baseline.EnsembleNaive
+
+    References:
+    ----------
+    [1]. https://otexts.com/fpp2/simple-methods.html
+
+    Examples:
+    --------
 
+    Basic fitting and prediction
+
+    >>> y_train = np.arange(10)
+    >>> model = Naive1()
+    >>> model.fit(y_train)
+    >>> model.predict(horizon=7)
+    array([9., 9., 9., 9., 9., 9., 9.]
+
+    fit_predict() convenience method
+
+    >>> y_train = np.arange(10)
+    >>> model = Naive1()
+    >>> model.fit_predict(y_train, horizon=7)
+    array([9., 9., 9., 9., 9., 9., 9.]
+
+    80 and 95% prediction intervals
+
+    >>> y_train = np.arange(10)
+    >>> model = Naive1()
+    >>> model.fit(y_train)
+    >>> y_pred, y_intervals = model.predict(horizon=2,
+                                            return_pred_interval=True,
+                                            alpha=[0.1, 0.05])
+    >>> y_pred
+    array([9., 9.]
+    >>> y_intervals[0]
+    array([[ 7.71844843, 10.28155157],
+           [ 7.1876124 , 10.8123876 ]])
+    >>> y_intervals[1]
+    array([[ 7.35514637, 10.64485363],
+           [ 6.67382569, 11.32617431]])
+
+
+    Fitted values (one step in-sample predictions)
+    .fittedvalue returns a pandas.Series called pred
+
+    >>> y_train = np.arange(5)
+    >>> model = Naive1()
+    >>> model.fit(y_train)
+    >>> model.fittedvalues
+    0    NaN
+    1    0.0
+    2    1.0
+    3    2.0
+    4    3.0
+    Name: pred, dtype: float64
+
+    '''
     def __init__(self):
         '''
         Constructor method
@@ -239,6 +315,18 @@ def __init__(self):
         '''
         self._fitted = None
 
+    def __repr__(self):
+        '''
+        String representation of object
+        '''
+        return f'Naive1()'
+
+    def __str__(self):
+        '''
+        Print/str representation of object
+        '''
+        return f'Naive1()'
+
     def fit(self, train):
         '''
         Train the naive model
@@ -305,7 +393,7 @@ def predict(self, horizon, return_predict_int=False, alpha=None):
             default behaviour is to return 80 and 90% intervals.
 
         Returns:
-        ------
+        -------
 
         if return_predict_int = False
 
@@ -352,9 +440,34 @@ class SNaive(Forecast):
     Each forecast to be equal to the last observed value from the
     same season of the year (e.g., the same month of the previous year).
 
-    SNF is useful for highly seasonal data.
+    SNF is useful for highly seasonal data. See [1]
 
-    See Hyndman: https://otexts.com/fpp2/simple-methods.html
+    Attributes
+    ----------
+    fittedvalues: pd.Series
+        In-sample predictions of training data
+    resid: pd.Series
+        In-sample residuals
+
+    Methods
+    -------
+    fit(train)
+        fit the model to training data
+    predict(horizon, return_predict_int=False, alpha=None)
+        Predict h-steps ahead
+    fit_predict(train, horizons, return_predict_int=False, alpha=None)
+        convenience method.  combine fit() and predict()
+
+    See Also
+    --------
+    forecast_tools.baseline.Naive1
+    forecast_tools.baseline.Drift
+    forecast_tools.baseline.Average
+    forecast_tools.baseline.EnsembleNaive
+
+    References:
+    -----------
+    [1]. https://otexts.com/fpp2/simple-methods.html
 
     '''
 
@@ -368,6 +481,18 @@ def __init__(self, period):
         self._period = period
         self._fitted = None
 
+    def __repr__(self):
+        '''
+        String representation of object
+        '''
+        return f'SNaive1(period={self._period})'
+
+    def __str__(self):
+        '''
+        Print/str representation of object
+        '''
+        return f'SNaive1(period={self._period})'
+
     def fit(self, train):
         '''
         Seasonal naive forecast - train the model
@@ -471,14 +596,52 @@ class Average(Forecast):
     Average forecast.  Forecast is set to the average
     of the historical data.
 
-    See Makridakis, Wheelwright and Hyndman (1998)
+    See for discussion of the average as a forecat measure [1]
+
+    Attributes
+    ----------
+    fittedvalues: pd.Series
+        In-sample predictions of training data
+    resid: pd.Series
+        In-sample residuals
+
+    Methods
+    -------
+    fit(train)
+        fit the model to training data
+    predict(horizon, return_predict_int=False, alpha=None)
+        Predict h-steps ahead
+    fit_predict(train, horizons, return_predict_int=False, alpha=None)
+        convenience method.  combine fit() and predict()
+
+    See Also
+    --------
+    forecast_tools.baseline.Naive1
+    forecast_tools.baseline.SNaive
+    forecast_tools.baseline.Drift
+    forecast_tools.baseline.EnsembleNaive
 
+    References:
+    -----------
+    [1.] Makridakis, Wheelwright and Hyndman. Forecasting (1998)
     '''
 
     def __init__(self):
         self._pred = None
         self._fitted = None
 
+    def __repr__(self):
+        '''
+        String representation of object
+        '''
+        return f'Average()'
+
+    def __str__(self):
+        '''
+        Print/str representation of object
+        '''
+        return f'Average()'
+
     def _get_fitted(self):
         return self._fitted['pred']
 
@@ -515,7 +678,6 @@ def fit(self, train):
         self._pred = _train.mean()
         # ddof set to get sample mean
         self._resid_std = (_train - self._pred).std(ddof=1)
-        print(self._resid_std)
         self._fitted['pred'] = self._pred
         self._fitted['resid'] = self._fitted['actual'] - self._fitted['pred']
 
@@ -582,23 +744,63 @@ def _std_h(self, horizon):
 
 class Drift(Forecast):
     '''
-    Naive1 with drift: Carry the last value foreward across a
-    forecast horizon but allow for upwards of downwards drift.
+    Naive1 forecast with drift
 
-    Drift = average change in the historical data.
+    Carry the last value foreward across a forecast horizon but
+    allow for upwards of downwards drift defined in [1]
 
-    https://otexts.com/fpp2/simple-methods.html
+    Drift = average change in the historical data.
 
     Note.  The current implementation has a standard error of the forecast
     that is the same as for the naive1 se.  This could be adjusted for drift.
     The following link suggests this is minor and benchmark with R is v.similar
-    https://www.coursehero.com/file/p12k3ln/For-the-random-walk-with-drift-model-the-1-step-ahead-forecast-standard-error/
+    [2]
+
+    Attributes
+    ----------
+    fittedvalues: pd.Series
+        In-sample predictions of training data
+    resid: pd.Series
+        In-sample residuals
+
+    Methods
+    -------
+    fit(train)
+        fit the model to training data
+    predict(horizon, return_predict_int=False, alpha=None)
+        Predict h-steps ahead
+    fit_predict(train, horizons, return_predict_int=False, alpha=None)
+        convenience method.  combine fit() and predict()
+
+    See Also
+    --------
+    forecast_tools.baseline.Naive1
+    forecast_tools.baseline.SNaive
+    forecast_tools.baseline.Average
+    forecast_tools.baseline.EnsembleNaive
+
+    References:
+    -----------
+    [1]. https://otexts.com/fpp2/simple-methods.html
+    [2]. https://www.coursehero.com/file/p12k3ln/For-the-random-walk-with-drift-model-the-1-step-ahead-forecast-standard-error/
 
     '''
 
     def __init__(self):
         self._fitted = None
 
+    def __repr__(self):
+        '''
+        String representation of object
+        '''
+        return f'Drift()'
+
+    def __str__(self):
+        '''
+        Print/str representation of object
+        '''
+        return f'Drift()'
+
     def _get_fitted_gradient(self):
         return self._fitted['gradient_fit']
 
@@ -702,7 +904,34 @@ def _std_h(self, horizon):
     fitted_gradient = property(_get_fitted_gradient)
 
 
-class EnsembleNaive(object):
+class EnsembleNaive(Forecast):
+    '''
+    An ensemble of all naive forecast methods.
+
+    Attributes
+    ----------
+    fittedvalues: pd.Series
+        In-sample predictions of training data
+    resid: pd.Series
+        In-sample residuals
+
+    Methods
+    -------
+    fit(train)
+        fit the model to training data
+    predict(horizon, return_predict_int=False, alpha=None)
+        Predict h-steps ahead
+    fit_predict(train, horizons, return_predict_int=False, alpha=None)
+        convenience method.  combine fit() and predict()
+
+    See Also
+    --------
+    forecast_tools.baseline.Naive1
+    forecast_tools.baseline.SNaive
+    forecast_tools.baseline.Average
+    forecast_tools.baseline.Drift
+    '''
+
     def __init__(self, seasonal_period):
         self._estimators = {'NF1': Naive1(),
                             'SNaive': SNaive(period=seasonal_period),
@@ -710,6 +939,20 @@ def __init__(self, seasonal_period):
                             'Drift': Drift()
                             }
 
+    def __repr__(self):
+        '''
+        String representation of object
+        '''
+        p = self._estimators['SNaive']._period
+        return f'EnsembleNaive(seasonal_period={p})'
+
+    def __str__(self):
+        '''
+        Print/str representation of object
+        '''
+        p = self._estimators['SNaive']._period
+        return f'EnsembleNaive(seasonal_period={p})'
+
     def fit(self, train):
         '''
         Parameters:
@@ -724,7 +967,7 @@ def fit(self, train):
         for _, estimator in self._estimators.items():
             estimator.fit(train)
 
-    def predict(self, horizon):
+    def predict(self, horizon, return_predict_int=False, alpha=None):
         preds = []
         for _, estimator in self._estimators.items():
             model_preds = estimator.predict(horizon)
@@ -732,6 +975,13 @@ def predict(self, horizon):
 
         return np.array(preds).mean(axis=0)
 
+    def _std_h(self, horizon):
+        '''
+        Calculate the standard error of the residuals over
+        a forecast horizon.  This is method specific.
+        '''
+        pass
+
 
 def baseline_estimators(seasonal_period):
     '''
diff --git a/forecast_tools/metrics.py b/forecast_tools/metrics.py
index 18ea642..092da88 100644
--- a/forecast_tools/metrics.py
+++ b/forecast_tools/metrics.py
@@ -7,12 +7,13 @@
 RMSE - root mean squared error
 MAPE - mean absolute percentage error
 sMAPE - symmetric MAPE.
+MASE - mean absolute scaled error
 
 coverage - prediction interval coverage
 '''
 import numpy as np
 
-from forecast_tools.baseline import Naive1, SNaive
+from forecast_tools.baseline import SNaive
 
 
 def as_arrays(y_true, y_pred):
@@ -179,6 +180,8 @@ def symmetric_mean_absolute_percentage_error(y_true, y_pred):
 
 def coverage(y_true, pred_intervals):
     '''
+    Prediction Interval Coverage
+
     Calculates the proportion of the true
     values are that are covered by the lower
     and upper bounds of the prediction intervals
@@ -231,18 +234,17 @@ def mean_absolute_scaled_error(y_true, y_pred, y_train, period=None):
     float,
         scalar value representing the MASE
     '''
+    y_true, y_pred = as_arrays(y_true, y_pred)
 
     if period is None:
-        in_sample = Naive1()
-        in_sample.fit(y_train)
+        period = 1
+
+    in_sample = SNaive(period=period)
+    in_sample.fit(y_train)
 
-    else:
-        in_sample = SNaive(period=period)
-        in_sample.fit(y_train)
-        y_train = y_train[period:]
+    mae_insample = mean_absolute_error(y_train[period:],
+                                       in_sample.fittedvalues.dropna())
 
-    mae_insample = mean_absolute_error(
-        y_train, in_sample.fittedvalues.dropna())
     return mean_absolute_error(y_true, y_pred) / mae_insample
 
 
@@ -271,14 +273,14 @@ def forecast_errors(y_true, y_pred, metrics='all'):
 
     Example:
     ---------
-    y_true = [45, 60, 23, 45]
-    y_preds = [50, 50, 50, 50]
+    >>> y_true = [45, 60, 23, 45]
+    >>> y_preds = [50, 50, 50, 50]
 
-    metrics = forecast_errors(y_true, y_preds)
-    print(metrics)
+    >>> metrics = forecast_errors(y_true, y_preds)
+    >>> print(metrics)
 
-    metrics = forecast_errors(y_true, y_preds, metrics=['mape', 'smape'])
-    print(metrics)
+    >>> metrics = forecast_errors(y_true, y_preds, metrics=['mape', 'smape'])
+    >>> print(metrics)
 
     '''
     y_true, y_pred = as_arrays(y_true, y_pred)
diff --git a/forecast_tools/model_selection.py b/forecast_tools/model_selection.py
index 62298a6..a6fc086 100644
--- a/forecast_tools/model_selection.py
+++ b/forecast_tools/model_selection.py
@@ -30,6 +30,10 @@
 import numpy as np
 from joblib import Parallel, delayed
 
+from forecast_tools.metrics import (mean_absolute_scaled_error,
+                                    _forecast_error_functions)
+from forecast_tools.baseline import baseline_estimators
+
 
 def rolling_forecast_origin(train, min_train_size, horizon, step=1):
     '''
@@ -95,6 +99,8 @@ def sliding_window(train, window_size, horizon, step=1):
 
 def forecast_accuracy(model, train, test, horizons, metric):
     '''
+    Forecast accuracy of a model over multiple horizons
+
     Returns forecast accuracy of model fit on
     training data and compared against a test set
     with a given metric.
@@ -134,19 +140,16 @@ def forecast_accuracy(model, train, test, horizons, metric):
     return h_accuracy
 
 
-def cross_validation_score(model, train, cv, metric, horizons=None,
+def cross_validation_score(model, cv, metric, horizons=None,
                            n_jobs=-1):
     '''
-    Calculate cross validation scores
+    Cross validation scores
 
     Parameters:
     ----------
     model: object
         forecast model
 
-    train: array-like
-        training data
-
     metric: func(y_true, y_pred)
         forecast error metric
 
@@ -183,3 +186,320 @@ def cross_validation_score(model, train, cv, metric, horizons=None,
                                     for cv_train, cv_test in cv)
 
     return np.array(cv_scores)
+
+
+def forecast(model, train, test, horizon):
+    '''
+    h-step prediction of a model
+
+    Returns a tuple of (y_preds, y_train, y_true) of model fit
+    to training data
+
+    Parameters:
+    ----------
+    model - object
+        forecasting model with .fit(train) and
+        .predict(horizon) methods
+
+    train - array-like
+        training data
+
+    test: array-like
+        holdout data for testing
+
+    horizon: int
+        forecast horizon
+
+    Returns:
+    --------
+    tuple (y_pred, y_train, y_true)
+    '''
+    y_pred = model.fit_predict(train, horizon)
+    return train, test, y_pred
+
+
+def cross_validation_folds(model, cv, n_jobs=-1):
+    '''
+    Cross validation forecasts
+
+    Parameters:
+    ----------
+    model: object
+        forecast model
+
+    cv: object
+        cross validation generator
+        i.e. rolling_forecast_origin or sliding_window
+
+    n_jobs: int, optional (default=-1)
+        when -1 runs across all cores
+        set = 1 to run each cross validation seperately.
+        using -1 speeds up cross validation of slow running models.
+
+    Returns:
+    -------
+    np.ndarray of tuples
+    each tuple is (cv_train, cv_test, cv_y_pred)
+    '''
+
+    cv_folds = \
+        Parallel(n_jobs=n_jobs)(delayed(forecast)(model,
+                                                  cv_train,
+                                                  cv_test,
+                                                  len(cv_test))
+                                for cv_train, cv_test in cv)
+
+    return np.array(cv_folds)
+
+
+def scaled_cross_validation_score(model, cv, seasonal_period=None):
+    '''
+    Mean absolute scaled error cross validation score
+
+    Parameters:
+    ----------
+    model: object
+        forecast model
+
+    cv: generator
+        time series cross validation fold generator
+
+    metric: func(y_true, y_py_red)
+        forecast error metric
+
+    seasonal_period: None or int, optional (default=None)
+        if none the in-sample one step Naive1 used for scaling.
+        if int SNaive is used instead.
+
+    Returns:
+    --------
+    np.array of mase scores
+    '''
+
+    folds = cross_validation_folds(model, cv)
+
+    scores = []
+    for y_train, y_true, y_pred in folds:
+        score = mean_absolute_scaled_error(y_true, y_pred, y_train,
+                                           period=seasonal_period)
+        scores.append(score)
+
+    return np.array(scores)
+
+
+def auto_naive(y_train, horizon=1, seasonal_period=1, min_train_size='auto',
+               method='cv', step=1, window_size='auto', metric='mae'):
+    '''Automatic selection of the 'best' naive benchmark on a 'single' series
+
+    The selection process uses out-of-sample cv performance.
+
+    By default auto_naive uses cross validation to estimate the mean
+    point forecast peformance of all naive methods.  It selects the method
+    with the lowest point forecast metric on average.
+
+    If there is limited data for training a basic holdout sample could be
+    used.
+
+    Dev note: the plan is to update this to work with multiple series.
+    It would be best to use MASE for multiple series comparison.
+
+    Parameters:
+    ----------
+    y_train: array-like
+        training data.  typically in a pandas.Series, pandas.DataFrame
+        or numpy.ndarray format.
+
+    horizon: int, optional (default=1)
+        Forecast horizon.
+
+    seasonal_period: int, optional (default=1)
+        Frequency of the data.  E.g. 7 for weekly pattern, 12 for monthly
+        365 for daily.
+
+    min_train_size: int or str, optional (default='auto')
+        The size of the initial training set (if method=='ro' or 'sw').
+        If 'auto' then then min_train_size is set to len(y_train) // 3
+        If main_train_size='auto' and method='holdout' then
+        min_train_size = len(y_train) - horizon.
+
+    method: str, optional (default='cv')
+        out of sample selection method.
+        'ro' - rolling forecast origin
+        'sw' - sliding window
+        'cv' - scores from both ro and sw
+        'holdout' - single train/test split
+         Methods'ro' and 'sw' are similar, however, sw has a fixed
+         window_size and drops older data from training.
+
+    step: int, optional (default=1)
+        The stride/step of the cross-validation. I.e. the number
+        of observations to move forward between folds.
+
+    window_size: str or int, optional (default='auto')
+        The window_size if using sliding window cross validation
+        When 'auto' and method='sw' then
+        window_size=len(y_train) // 3
+
+    metric: str, optional (default='mae')
+        The metric to measure out of sample accuracy.
+        Options: mase, mae, mape, smape, mse, rmse, me.
+
+    Returns:
+    --------
+    dict
+        'model': baseline.Forecast
+        f'{metric}': float
+
+        Contains the model and its CV performance.
+
+    Raises:
+    -------
+    ValueError
+        For invalid method, metric, window_size parameters
+
+    See Also:
+    --------
+    forecast_tools.baseline.Naive1
+    forecast_tools.baseline.SNaive
+    forecast_tools.baseline.Drift
+    forecast_tools.baseline.Average
+    forecast_tools.baseline.EnsembleNaive
+    forecast_tools.baseline.baseline_estimators
+    forecast_tools.model_selection.rolling_forecast_origin
+    forecast_tools.model_selection.sliding_window
+    forecast_tools.model_selection.mase_cross_validation_score
+    forecast_tools.metrics.mean_absolute_scaled_error
+
+    Examples:
+    ---------
+    Measuring MAE and taking the best method using both
+    rolling origin and sliding window cross validation
+    of a 56 day forecast.
+
+    >>> from forecast_tools.datasets import load_emergency_dept
+    >>> y_train = load_emergency_dept
+    >>> best = auto_naive(y_train, seasonal_period=7, horizon=56)
+    >>> best
+    {'model': Average(), 'mae': 19.63791579700355}
+
+
+    Take a step of 7 days between cv folds.
+
+    >>> from forecast_tools.datasets import load_emergency_dept
+    >>> y_train = load_emergency_dept
+    >>> best = auto_naive(y_train, seasonal_period=7, horizon=56,
+        ...               step=7)
+    >>> best
+    {'model': Average(), 'mae': 19.675635558539383}
+
+    '''
+    valid_methods = ['holdout', 'ro', 'sw', 'cv']
+    metrics = _forecast_error_functions()
+    metrics['mase'] = mean_absolute_scaled_error
+
+    if method not in valid_methods:
+        raise ValueError(f"Method must be in {valid_methods}")
+
+    if metric not in metrics:
+        raise ValueError(f"Please select a metric from {list(metrics.keys())}")
+
+    if min_train_size == 'auto':
+        min_train_size = len(y_train) // 3
+    elif not type(min_train_size) is int:
+        raise ValueError(f"valid min_train_size values are 'auto' or int > 0")
+    elif min_train_size < 1:
+        raise ValueError(f"valid min_train_size values are 'auto' or int > 0")
+
+    if window_size == 'auto':
+        window_size = len(y_train) // 3
+    elif not type(window_size) is int:
+        raise ValueError(f"valid window_size values are 'auto' or int > 0")
+    elif window_size < 1:
+        raise ValueError(f"valid window_size values are 'auto' or int > 0")
+
+    baselines = baseline_estimators(seasonal_period)
+
+    method_score = []
+    if method == 'cv':
+        for _, model in baselines.items():
+            cv_ro = rolling_forecast_origin(train=y_train,
+                                            min_train_size=min_train_size,
+                                            horizon=horizon,
+                                            step=step)
+
+            cv_sw = sliding_window(train=y_train,
+                                   window_size=window_size,
+                                   horizon=horizon,
+                                   step=step)
+
+            if metric == 'mase':
+                score_ro = scaled_cross_validation_score(model, cv_ro,
+                                                         seasonal_period)
+                score_sw = scaled_cross_validation_score(model, cv_sw,
+                                                         seasonal_period)
+
+            else:
+                score_ro = cross_validation_score(model, cv_ro,
+                                                  metrics[metric])
+                score_sw = cross_validation_score(model, cv_sw,
+                                                  metrics[metric])
+
+            score = np.concatenate([score_ro, score_sw])
+            method_score.append(score.mean())
+
+    elif method == 'ro':
+        for _, model in baselines.items():
+            cv = rolling_forecast_origin(train=y_train,
+                                         min_train_size=min_train_size,
+                                         horizon=horizon,
+                                         step=step)
+
+            if metric == 'mase':
+                score_ro = scaled_cross_validation_score(model, cv,
+                                                         seasonal_period)
+
+            else:
+                score_ro = cross_validation_score(model, cv, metrics[metric])
+
+            method_score.append(score_ro.mean())
+
+    elif method == 'sw':
+        for _, model in baselines.items():
+            cv = sliding_window(train=y_train,
+                                window_size=window_size,
+                                horizon=horizon,
+                                step=step)
+
+            if metric == 'mase':
+                score_sw = scaled_cross_validation_score(model, cv,
+                                                         seasonal_period)
+
+            else:
+                score_sw = cross_validation_score(model, cv, metrics[metric])
+
+            method_score.append(score_sw.mean())
+
+    else:
+        # single train test split
+        min_train_size = len(y_train) - horizon
+        train = y_train[:min_train_size]
+        test = y_train[min_train_size:]
+
+        for _, model in baselines.items():
+            model.fit(train)
+            y_preds = model.predict(horizon)
+            if metric == 'mase':
+                score = metrics[metric](test, y_preds, y_train,
+                                        seasonal_period)
+            else:
+                score = metrics[metric](test, y_preds)
+
+            method_score.append(score.mean())
+
+    method_score = np.array(method_score)
+    best_index = np.argmin(method_score)
+
+    best = {'model': list(baselines.items())[best_index][1],
+            f'{metric}': method_score[best_index]}
+
+    return best
diff --git a/setup.py b/setup.py
index 7265eaa..df4de89 100644
--- a/setup.py
+++ b/setup.py
@@ -12,7 +12,7 @@
 
 setuptools.setup(
     name="forecast-tools",
-    version="0.1.4",
+    version="0.1.5",
     author="Thomas Monks",
     author_email="forecast_tools@gmail.com",
     license="The MIT License (MIT)",


Test patch:
diff --git a/tests/test_metrics.py b/tests/test_metrics.py
index 959e8f4..de0882c 100644
--- a/tests/test_metrics.py
+++ b/tests/test_metrics.py
@@ -169,3 +169,34 @@ def test_coverage(y_true, y_intervals, expected):
     y_intervals = np.array(y_intervals).T
     error = m.coverage(y_true, y_intervals)
     assert pytest.approx(expected) == error
+
+
+@pytest.mark.parametrize("y_train, y_pred, y_true, expected",
+                         [(np.arange(10), [1, 2, 3, 4, 5, 6],
+                          [1, 2, 3, 4, 5, 6], 0.0),
+                          (np.arange(1, 21), np.arange(21, 26),
+                           np.full(5, 10), 13)])
+def test_mase_naive(y_train, y_true, y_pred, expected):
+    '''
+    test mean absolute scaled error calculation using naive as scaler.
+
+    test calcs produced using libre office calc.
+    '''
+    error = m.mean_absolute_scaled_error(y_true, y_pred, y_train, period=None)
+    assert pytest.approx(expected) == error
+
+
+@pytest.mark.parametrize("y_train, y_pred, y_true, expected",
+                         [(np.arange(1, 21), [1, 2, 3, 4, 5, 6],
+                          [1, 2, 3, 4, 5, 6], 0.0),
+
+                          (np.arange(1, 21), np.arange(21, 26),
+                           np.full(5, 10), 1.85714286)])
+def test_mase_snaive(y_train, y_true, y_pred, expected):
+    '''
+    test mean absolute scaled error calculation using SNaive as scaler
+
+    test calcs produced using libre office calc.
+    '''
+    error = m.mean_absolute_scaled_error(y_true, y_pred, y_train, period=7)
+    assert pytest.approx(expected) == error
diff --git a/tests/test_model_selection.py b/tests/test_model_selection.py
index 2d8c08b..6f85bf8 100644
--- a/tests/test_model_selection.py
+++ b/tests/test_model_selection.py
@@ -248,3 +248,116 @@ def test_rfo_number_of_folds_pd(train_size, min_train_size, horizon, step,
         actual += 1
 
     assert expected == actual
+
+
+@pytest.mark.parametrize("train_size, min_train_size, horizon, step, expected",
+                         [(10, 3, 1, 1, 7),
+                          (10, 3, 2, 1, 6),
+                          (10, 3, 1, 2, 4),
+                          (10, 3, 2, 2, 3)])
+def test_mase_cv_number_of_folds(train_size, min_train_size, horizon, step,
+                                 expected):
+    '''
+    check that the number of folds returned from rolling origin
+    is as expected when data source is a pandas.DataFrame
+    '''
+    train = np.arange(train_size)
+    cv = ms.rolling_forecast_origin(train, min_train_size=min_train_size,
+                                    horizon=horizon, step=step)
+
+    scores = ms.scaled_cross_validation_score(b.Naive1(), cv)
+
+    print(expected, len(scores))
+    print(scores)
+    assert expected == len(scores)
+
+
+@pytest.mark.parametrize("train, metric",
+                         [(np.arange(7), 'mae'),
+                          (np.arange(7), 'mase')])
+def test_auto_naive_return_length(train, metric):
+    result = ms.auto_naive(train, metric=metric)
+    assert len(result) == 2
+
+
+@pytest.mark.parametrize("train, cv, metric",
+                         [(np.arange(7), 'ro', 'mase'),
+                          (np.arange(7), 'sw', 'mase'),
+                          (np.arange(7), 'holdout', 'mase'),
+                          (np.arange(7), 'ro', 'mae'),
+                          (np.arange(7), 'sw', 'mae'),
+                          (np.arange(7), 'holdout', 'mae')])
+def test_auto_naive_return_length_vary_cv(train, cv, metric):
+    result = ms.auto_naive(train, method=cv, metric=metric)
+    assert len(result) == 2
+
+
+@pytest.mark.parametrize("metric",
+                         [('ma'),
+                          (999),
+                          ('xxx')])
+def test_auto_naive_metric_value_error(metric):
+    train = np.arange(5)
+    with pytest.raises(ValueError):
+        ms.auto_naive(train, metric=metric)
+
+
+@pytest.mark.parametrize("method",
+                         [('ma'),
+                          (999),
+                          ('xxx')])
+def test_auto_naive_method_value_error(method):
+    train = np.arange(5)
+    with pytest.raises(ValueError):
+        ms.auto_naive(train, method=method)
+
+
+@pytest.mark.parametrize("window_size",
+                         [('ma'),
+                          (-999),
+                          ('xxx')])
+def test_auto_naive_window_size_value_error(window_size):
+    train = np.arange(5)
+    with pytest.raises(ValueError):
+        ms.auto_naive(train, window_size=window_size)
+
+
+@pytest.mark.parametrize("min_train_size",
+                         [('ma'),
+                          (-999),
+                          ('xxx')])
+def test_auto_naive_mintrain_size_value_error(min_train_size):
+    train = np.arange(5)
+    with pytest.raises(ValueError):
+        ms.auto_naive(train, min_train_size=min_train_size)
+
+
+@pytest.mark.parametrize("train_size, min_train_size, horizon, step, expected",
+                         [(10, 3, 1, 1, 7),
+                          (10, 3, 2, 1, 6),
+                          (10, 3, 1, 2, 4),
+                          (10, 3, 2, 2, 3)])
+def test_cross_validation_folds_return_length(train_size, min_train_size,
+                                              horizon, step, expected):
+    y_train = np.arange(train_size)
+    cv = ms.rolling_forecast_origin(y_train, min_train_size, horizon, step)
+    result = ms.cross_validation_folds(b.Naive1(), cv)
+
+    assert len(result) == expected
+
+
+@pytest.mark.parametrize("train_size, min_train_size, horizon, step, expected",
+                         [(10, 3, 1, 1, 7),
+                          (10, 3, 2, 1, 6),
+                          (10, 3, 1, 2, 4),
+                          (10, 3, 2, 2, 3)])
+def test_cross_validation_score_return_length(train_size, min_train_size,
+                                              horizon, step, expected):
+    y_train = np.arange(train_size)
+    cv = ms.rolling_forecast_origin(y_train, min_train_size, horizon, step)
+    metric = metrics.mean_error
+    result = ms.cross_validation_score(b.Naive1(), cv, metric=metric, n_jobs=1)
+
+    assert len(result) == expected
+
+