enh(SES): improve implementation

README.md

@@ -170,10 +170,10 @@ Uses a weighted average of all past observations where the weights decrease expo
 
 |Model | Point Forecast | Probabilistic Forecast | Insample fitted values | Probabilistic fitted values |Exogenous features|
 |:------|:-------------:|:----------------------:|:---------------------:|:----------------------------:|:----------------:|
-|[SimpleExponentialSmoothing](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#simpleexponentialsmoothing)|✅|||||
-|[SimpleExponentialSmoothingOptimized](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#simpleexponentialsmoothingoptimized)|✅|||||
-|[SeasonalExponentialSmoothing](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#seasonalexponentialsmoothing)|✅|||||
-|[SeasonalExponentialSmoothingOptimized](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#seasonalexponentialsmoothingoptimized)|✅|||||
+|[SimpleExponentialSmoothing](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#simpleexponentialsmoothing)|✅||✅|||
+|[SimpleExponentialSmoothingOptimized](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#simpleexponentialsmoothingoptimized)|✅||✅|||
+|[SeasonalExponentialSmoothing](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#seasonalexponentialsmoothing)|✅||✅|||
+|[SeasonalExponentialSmoothingOptimized](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#seasonalexponentialsmoothingoptimized)|✅||✅|||
 |[Holt](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#holt)|✅|✅|✅|✅||
 |[HoltWinters](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#holtwinters)|✅|✅|✅|✅||
 

nbs/src/core/models.ipynb

@@ -58,11 +58,11 @@
     "#| export\n",
     "import warnings\n",
     "from math import trunc\n",
-    "from typing import Any, Dict, List, Optional, Sequence, Tuple, Union\n",
+    "from typing import Any, Dict, List, Optional, Tuple, Union\n",
     "\n",
     "import numpy as np\n",
     "from numba import njit\n",
-    "from scipy.optimize import minimize\n",
+    "from scipy.optimize import minimize_scalar\n",
     "from scipy.special import inv_boxcox\n",
     "\n",
     "from statsforecast.arima import (\n",
@@ -3649,37 +3649,61 @@
    "source": [
     "#| export\n",
     "@njit(nogil=NOGIL, cache=CACHE)\n",
-    "def _ses_fcst_mse(x: np.ndarray, alpha: float) -> Tuple[float, float, np.ndarray]:\n",
-    "    r\"\"\"Perform simple exponential smoothing on a series.\n",
+    "def _ses_sse(alpha: float, x: np.ndarray) -> float:\n",
+    "    r\"\"\"Compute the residual sum of squares for a simple exponential smoothing fit.\n",
+    "    \n",
+    "    Parameters\n",
+    "    ----------\n",
+    "    alpha : float \n",
+    "        Smoothing parameter.\n",
+    "    x : numpy.array\n",
+    "        Clean time series of shape (n, ).\n",
+    "\n",
+    "    Returns\n",
+    "    -------\n",
+    "    sse : float\n",
+    "        Residual sum of squares for the fit.\n",
     "\n",
-    "    This function returns the one step ahead prediction\n",
-    "    as well as the mean squared error of the fit.\n",
     "    \"\"\"\n",
-    "    smoothed = x[0]\n",
-    "    n = x.size\n",
-    "    mse = 0.\n",
-    "    fitted = np.full(n, np.nan, dtype=x.dtype)\n",
+    "    complement = 1 - alpha\n",
+    "    forecast = x[0]\n",
+    "    sse = 0.0\n",
+    "    \n",
+    "    for i in range(1, len(x)):\n",
+    "        forecast = alpha * x[i - 1] + complement * forecast\n",
+    "        sse += (x[i] - forecast) ** 2\n",
+    "    \n",
+    "    return sse\n",
     "\n",
-    "    for i in range(1, n):\n",
-    "        smoothed = (alpha * x[i - 1] + (1 - alpha) * smoothed).item()\n",
-    "        error = x[i] - smoothed\n",
-    "        mse += error * error\n",
-    "        fitted[i] = smoothed\n",
     "\n",
-    "    mse /= n\n",
-    "    forecast = alpha * x[-1] + (1 - alpha) * smoothed\n",
-    "    return forecast, mse, fitted\n",
+    "@njit(nogil=NOGIL, cache=CACHE)\n",
+    "def _ses_forecast(x: np.ndarray, alpha: float) -> Tuple[float, np.ndarray]:\n",
+    "    r\"\"\"Compute the one-step ahead forecast for a simple exponential smoothing fit.\n",
     "\n",
+    "    Parameters\n",
+    "    ----------\n",
+    "    x : numpy.array\n",
+    "        Clean time series of shape (n, ).\n",
+    "    alpha : float\n",
+    "        Smoothing parameter.\n",
     "\n",
-    "def _ses_mse(alpha: float, x: np.ndarray) -> float:\n",
-    "    r\"\"\"Compute the mean squared error of a simple exponential smoothing fit.\"\"\"\n",
-    "    _, mse, _ = _ses_fcst_mse(x, alpha)\n",
-    "    return mse\n",
+    "    Returns\n",
+    "    -------\n",
+    "    tuple of (float, numpy.array)\n",
+    "        One-step ahead forecast and in-sample fitted values.\n",
     "\n",
+    "    \"\"\"\n",
+    "    complement = 1 - alpha\n",
+    "    fitted = np.empty_like(x)\n",
+    "    fitted[0] = x[0]\n",
+    "    j = 0\n",
     "\n",
-    "def _ses_forecast(x: np.ndarray, alpha: float) -> Tuple[float, np.ndarray]:\n",
-    "    r\"\"\"One step ahead forecast with simple exponential smoothing.\"\"\"\n",
-    "    forecast, _, fitted = _ses_fcst_mse(x, alpha)\n",
+    "    for i in range(1, len(x)):\n",
+    "        fitted[i] = alpha * x[j] + complement * fitted[j]\n",
+    "        j += 1\n",
+    "\n",
+    "    forecast = alpha * x[j] + complement * fitted[j]\n",
+    "    fitted[0] = np.nan\n",
     "    return forecast, fitted\n",
     "\n",
     "\n",
@@ -3701,16 +3725,28 @@
     "\n",
     "def _optimized_ses_forecast(\n",
     "        x: np.ndarray,\n",
-    "        bounds: Sequence[Tuple[float, float]] = [(0.1, 0.3)]\n",
+    "        bounds: Tuple[float, float] = (0.1, 0.3)\n",
     "    ) -> Tuple[float, np.ndarray]:\n",
-    "    r\"\"\"Searches for the optimal alpha and computes SES one step forecast.\"\"\"\n",
-    "    alpha = minimize(\n",
-    "        fun=_ses_mse,\n",
-    "        x0=(0,),\n",
-    "        args=(x,),\n",
+    "    r\"\"\"Compute the one-step ahead forecast for an optimal simple exponential smoothing fit.\n",
+    "    \n",
+    "    Parameters\n",
+    "    ----------\n",
+    "    x : numpy.array\n",
+    "        Clean time series of shape (n, ).\n",
+    "    bounds : tuple of (float, float)\n",
+    "        Lower and upper optimisation bounds for alpha.\n",
+    "\n",
+    "    Returns\n",
+    "    -------\n",
+    "    tuple of (float, numpy.array)\n",
+    "        One-step ahead forecast and in-sample fitted values.\n",
+    "\n",
+    "    \"\"\"\n",
+    "    alpha = minimize_scalar(\n",
+    "        fun=_ses_sse,\n",
     "        bounds=bounds,\n",
-    "        method='L-BFGS-B'\n",
-    "    ).x[0]\n",
+    "        args=(x,),\n",
+    "    ).x\n",
     "    forecast, fitted = _ses_forecast(x, alpha)\n",
     "    return forecast, fitted\n",
     "\n",
@@ -3737,11 +3773,11 @@
     "    fitted: bool, # fitted values\n",
     "    alpha: float, # smoothing parameter\n",
     ") -> Dict[str, np.ndarray]: \n",
-    "    fcst, _, fitted_vals = _ses_fcst_mse(y, alpha)\n",
-    "    fcst = {'mean': _repeat_val(val=fcst, h=h)}\n",
+    "    fcst, fitted_vals = _ses_forecast(y, alpha)\n",
+    "    out = {'mean': _repeat_val(val=fcst, h=h)}\n",
     "    if fitted:\n",
-    "        fcst['fitted'] = fitted_vals\n",
-    "    return fcst"
+    "        out['fitted'] = fitted_vals\n",
+    "    return out"
    ]
   },
   {
@@ -4057,7 +4093,7 @@
     "        h: int, # forecasting horizon\n",
     "        fitted: bool, # fitted values\n",
     "    ):\n",
-    "    fcst_, fitted_vals = _optimized_ses_forecast(y, [(0.01, 0.99)])\n",
+    "    fcst_, fitted_vals = _optimized_ses_forecast(y, (0.01, 0.99))\n",
     "    mean = _repeat_val(val=fcst_, h=h)\n",
     "    fcst = {'mean': mean}\n",
     "    if fitted:\n",
@@ -4746,7 +4782,7 @@
     "    fitted_vals = np.full_like(y, np.nan)\n",
     "    for i in range(season_length):\n",
     "        init_idx = (i + n % season_length)\n",
-    "        season_vals[i], fitted_vals[init_idx::season_length] = _optimized_ses_forecast(y[init_idx::season_length], [(0.01, 0.99)])\n",
+    "        season_vals[i], fitted_vals[init_idx::season_length] = _optimized_ses_forecast(y[init_idx::season_length], (0.01, 0.99))\n",
     "    out = _repeat_val_seas(season_vals=season_vals, h=h)\n",
     "    fcst = {'mean': out}\n",
     "    if fitted:\n",
@@ -13027,9 +13063,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "python3",
+   "display_name": "statsforecast",
    "language": "python",
-   "name": "python3"
+   "name": "statsforecast"
   }
  },
  "nbformat": 4,

python/statsforecast/_modidx.py

@@ -720,13 +720,11 @@
                                       'statsforecast.models._seasonal_window_average': ( 'src/core/models.html#_seasonal_window_average',
                                                                                          'statsforecast/models.py'),
                                       'statsforecast.models._ses': ('src/core/models.html#_ses', 'statsforecast/models.py'),
-                                      'statsforecast.models._ses_fcst_mse': ( 'src/core/models.html#_ses_fcst_mse',
-                                                                              'statsforecast/models.py'),
                                       'statsforecast.models._ses_forecast': ( 'src/core/models.html#_ses_forecast',
                                                                               'statsforecast/models.py'),
-                                      'statsforecast.models._ses_mse': ('src/core/models.html#_ses_mse', 'statsforecast/models.py'),
                                       'statsforecast.models._ses_optimized': ( 'src/core/models.html#_ses_optimized',
                                                                                'statsforecast/models.py'),
+                                      'statsforecast.models._ses_sse': ('src/core/models.html#_ses_sse', 'statsforecast/models.py'),
                                       'statsforecast.models._tsb': ('src/core/models.html#_tsb', 'statsforecast/models.py'),
                                       'statsforecast.models._window_average': ( 'src/core/models.html#_window_average',
                                                                                 'statsforecast/models.py')},

python/statsforecast/models.py

@@ -12,11 +12,11 @@
 # %% ../../nbs/src/core/models.ipynb 4
 import warnings
 from math import trunc
-from typing import Any, Dict, List, Optional, Sequence, Tuple, Union
+from typing import Any, Dict, List, Optional, Tuple, Union
 
 import numpy as np
 from numba import njit
-from scipy.optimize import minimize
+from scipy.optimize import minimize_scalar
 from scipy.special import inv_boxcox
 
 from statsforecast.arima import (
@@ -2130,37 +2130,61 @@ def __init__(
 
 # %% ../../nbs/src/core/models.ipynb 129
 @njit(nogil=NOGIL, cache=CACHE)
-def _ses_fcst_mse(x: np.ndarray, alpha: float) -> Tuple[float, float, np.ndarray]:
-    r"""Perform simple exponential smoothing on a series.
+def _ses_sse(alpha: float, x: np.ndarray) -> float:
+    r"""Compute the residual sum of squares for a simple exponential smoothing fit.
 
-    This function returns the one step ahead prediction
-    as well as the mean squared error of the fit.
-    """
-    smoothed = x[0]
-    n = x.size
-    mse = 0.0
-    fitted = np.full(n, np.nan, dtype=x.dtype)
+    Parameters
+    ----------
+    alpha : float
+        Smoothing parameter.
+    x : numpy.array
+        Clean time series of shape (n, ).
 
-    for i in range(1, n):
-        smoothed = (alpha * x[i - 1] + (1 - alpha) * smoothed).item()
-        error = x[i] - smoothed
-        mse += error * error
-        fitted[i] = smoothed
+    Returns
+    -------
+    sse : float
+        Residual sum of squares for the fit.
 
-    mse /= n
-    forecast = alpha * x[-1] + (1 - alpha) * smoothed
-    return forecast, mse, fitted
+    """
+    complement = 1 - alpha
+    forecast = x[0]
+    sse = 0.0
 
+    for i in range(1, len(x)):
+        forecast = alpha * x[i - 1] + complement * forecast
+        sse += (x[i] - forecast) ** 2
 
-def _ses_mse(alpha: float, x: np.ndarray) -> float:
-    r"""Compute the mean squared error of a simple exponential smoothing fit."""
-    _, mse, _ = _ses_fcst_mse(x, alpha)
-    return mse
+    return sse
 
 
+@njit(nogil=NOGIL, cache=CACHE)
 def _ses_forecast(x: np.ndarray, alpha: float) -> Tuple[float, np.ndarray]:
-    r"""One step ahead forecast with simple exponential smoothing."""
-    forecast, _, fitted = _ses_fcst_mse(x, alpha)
+    r"""Compute the one-step ahead forecast for a simple exponential smoothing fit.
+
+    Parameters
+    ----------
+    x : numpy.array
+        Clean time series of shape (n, ).
+    alpha : float
+        Smoothing parameter.
+
+    Returns
+    -------
+    tuple of (float, numpy.array)
+        One-step ahead forecast and in-sample fitted values.
+
+    """
+    complement = 1 - alpha
+    fitted = np.empty_like(x)
+    fitted[0] = x[0]
+    j = 0
+
+    for i in range(1, len(x)):
+        fitted[i] = alpha * x[j] + complement * fitted[j]
+        j += 1
+
+    forecast = alpha * x[j] + complement * fitted[j]
+    fitted[0] = np.nan
     return forecast, fitted
 
 
@@ -2181,12 +2205,28 @@ def _probability(x: np.ndarray) -> np.ndarray:
 
 
 def _optimized_ses_forecast(
-    x: np.ndarray, bounds: Sequence[Tuple[float, float]] = [(0.1, 0.3)]
+    x: np.ndarray, bounds: Tuple[float, float] = (0.1, 0.3)
 ) -> Tuple[float, np.ndarray]:
-    r"""Searches for the optimal alpha and computes SES one step forecast."""
-    alpha = minimize(
-        fun=_ses_mse, x0=(0,), args=(x,), bounds=bounds, method="L-BFGS-B"
-    ).x[0]
+    r"""Compute the one-step ahead forecast for an optimal simple exponential smoothing fit.
+
+    Parameters
+    ----------
+    x : numpy.array
+        Clean time series of shape (n, ).
+    bounds : tuple of (float, float)
+        Lower and upper optimisation bounds for alpha.
+
+    Returns
+    -------
+    tuple of (float, numpy.array)
+        One-step ahead forecast and in-sample fitted values.
+
+    """
+    alpha = minimize_scalar(
+        fun=_ses_sse,
+        bounds=bounds,
+        args=(x,),
+    ).x
     forecast, fitted = _ses_forecast(x, alpha)
     return forecast, fitted
 
@@ -2206,11 +2246,11 @@ def _ses(
     fitted: bool,  # fitted values
     alpha: float,  # smoothing parameter
 ) -> Dict[str, np.ndarray]:
-    fcst, _, fitted_vals = _ses_fcst_mse(y, alpha)
-    fcst = {"mean": _repeat_val(val=fcst, h=h)}
+    fcst, fitted_vals = _ses_forecast(y, alpha)
+    out = {"mean": _repeat_val(val=fcst, h=h)}
     if fitted:
-        fcst["fitted"] = fitted_vals
-    return fcst
+        out["fitted"] = fitted_vals
+    return out
 
 # %% ../../nbs/src/core/models.ipynb 131
 class SimpleExponentialSmoothing(_TS):
@@ -2380,7 +2420,7 @@ def _ses_optimized(
     h: int,  # forecasting horizon
     fitted: bool,  # fitted values
 ):
-    fcst_, fitted_vals = _optimized_ses_forecast(y, [(0.01, 0.99)])
+    fcst_, fitted_vals = _optimized_ses_forecast(y, (0.01, 0.99))
     mean = _repeat_val(val=fcst_, h=h)
     fcst = {"mean": mean}
     if fitted:
@@ -2759,7 +2799,7 @@ def _seasonal_ses_optimized(
     for i in range(season_length):
         init_idx = i + n % season_length
         season_vals[i], fitted_vals[init_idx::season_length] = _optimized_ses_forecast(
-            y[init_idx::season_length], [(0.01, 0.99)]
+            y[init_idx::season_length], (0.01, 0.99)
         )
     out = _repeat_val_seas(season_vals=season_vals, h=h)
     fcst = {"mean": out}