add msarimax

README.md

@@ -110,6 +110,7 @@ Automatic forecasting tools search for the best parameters and select the best p
 |Model | Point Forecast | Probabilistic Forecast | Insample fitted values | Probabilistic fitted values |Exogenous features|
 |:------|:-------------:|:----------------------:|:---------------------:|:----------------------------:|:----------------:|
 |[AutoARIMA](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#autoarima)|✅|✅|✅|✅|✅|
+|[AutoMSARIMAX](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#automsarimax)|✅|✅|✅|✅|✅|
 |[AutoETS](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#autoets)|✅|✅|✅|✅||
 |[AutoCES](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#autoces)|✅|✅|✅|✅||
 |[AutoTheta](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#autotheta)|✅|✅|✅|✅||
@@ -123,6 +124,7 @@ These models exploit the existing autocorrelations in the time series.
 |Model | Point Forecast | Probabilistic Forecast | Insample fitted values | Probabilistic fitted values |Exogenous features|
 |:------|:-------------:|:----------------------:|:---------------------:|:----------------------------:|:----------------:|
 |[ARIMA](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#arima)|✅|✅|✅|✅|✅|
+|[MSARIMAX](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#msarimax)|✅|✅|✅|✅|✅|
 |[AutoRegressive](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#autoregressive)|✅|✅|✅|✅|✅|
 
 ### Theta Family
@@ -142,6 +144,7 @@ Suited for signals with more than one clear seasonality. Useful for low-frequenc
 
 |Model | Point Forecast | Probabilistic Forecast | Insample fitted values | Probabilistic fitted values |Exogenous features|
 |:------|:-------------:|:----------------------:|:---------------------:|:----------------------------:|:----------------:|
+|[MSARIMAX](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#msarimax)|✅|✅|✅|✅|✅|
 |[MSTL](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#mstl)|✅|✅|✅|✅|If trend forecaster supports|
 |[MFLES](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#mfles)|✅|✅|✅|✅|✅|
 |[TBATS](https://nixtlaverse.nixtla.io/statsforecast/src/core/models.html#tbats)|✅|✅|✅|✅||

docs/mintlify/docs.json

@@ -76,6 +76,7 @@
           "docs/models/arch.html",
           "docs/models/arima.html",
           "docs/models/autoarima.html",
+          "docs/models/automsarimax.html",
           "docs/models/autoces.html",
           "docs/models/autoets.html",
           "docs/models/autoregressive.html",
@@ -90,6 +91,7 @@
           "docs/models/holtwinters.html",
           "docs/models/imapa.html",
           "docs/models/mfles.html",
+          "docs/models/msarimax.html",
           "docs/models/multipleseasonaltrend.html",
           "docs/models/optimizedtheta.html",
           "docs/models/seasonalexponentialsmoothing.html",
@@ -133,4 +135,4 @@
       "href": "https://github.com/Nixtla/nixtla"
     }
   }
-}
+}

docs/src/core/models.html.md

@@ -7,7 +7,7 @@ title: Models
 
 StatsForecast offers a wide variety of statistical forecasting models grouped into the following categories:
 
-- **Auto Forecast**: Automatic forecasting tools that search for the best parameters and select the best possible model. Useful for large collections of univariate time series. Includes: AutoARIMA, AutoETS, AutoTheta, AutoCES, AutoMFLES, AutoTBATS.
+- **Auto Forecast**: Automatic forecasting tools that search for the best parameters and select the best possible model. Useful for large collections of univariate time series. Includes: AutoARIMA, AutoMSARIMAX, AutoETS, AutoTheta, AutoCES, AutoMFLES, AutoTBATS.
 
 - **ARIMA Family**: AutoRegressive Integrated Moving Average models for capturing autocorrelations in time series data.
 
@@ -17,7 +17,7 @@ StatsForecast offers a wide variety of statistical forecasting models grouped in
 
 - **Sparse or Intermittent**: Models suited for series with very few non-zero observations: ADIDA, CrostonClassic, CrostonOptimized, CrostonSBA, IMAPA, TSB.
 
-- **Multiple Seasonalities**: Models suited for signals with more than one clear seasonality. Useful for low-frequency data like electricity and logs: MSTL, MFLES, TBATS.
+- **Multiple Seasonalities**: Models suited for signals with more than one clear seasonality. Useful for low-frequency data like electricity and logs: MSARIMAX, MSTL, MFLES, TBATS.
 
 - **Theta Models**: Fit two theta lines to a deseasonalized time series using different techniques: Theta, OptimizedTheta, DynamicTheta, DynamicOptimizedTheta.
 
@@ -40,6 +40,19 @@ StatsForecast offers a wide variety of statistical forecasting models grouped in
         - predict_in_sample
         - forecast
 
+### AutoMSARIMAX
+
+::: statsforecast.models.AutoMSARIMAX
+    options:
+      show_source: true
+      heading_level: 4
+      members:
+        - __init__
+        - fit
+        - predict
+        - predict_in_sample
+        - forecast
+
 ### AutoETS
 
 ::: statsforecast.models.AutoETS
@@ -367,6 +380,19 @@ StatsForecast offers a wide variety of statistical forecasting models grouped in
 
 ## Multiple Seasonalities
 
+### MSARIMAX
+
+::: statsforecast.models.MSARIMAX
+    options:
+      show_source: true
+      heading_level: 4
+      members:
+        - __init__
+        - fit
+        - predict
+        - predict_in_sample
+        - forecast
+
 ### MSTL
 
 ::: statsforecast.models.MSTL

nbs/docs/models/AutoMSARIMAX.ipynb

@@ -0,0 +1,168 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "594abb0a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#| echo: false\n",
+    "import warnings\n",
+    "warnings.filterwarnings(\"ignore\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# AutoMSARIMAX Model\n",
+    "\n",
+    "> Automatic multiple-seasonal SARIMAX model.\n",
+    "\n",
+    "`AutoMSARIMAX` searches over bounded AR, differencing, and MA orders at each lag level, then selects the specification with the best information criterion. It fits one SARIMAX model for every candidate order combination, so broad bounds can be slow. Start with restricted bounds and widen them only when needed.\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Simulated Data\n",
+    "\n",
+    "The example below simulates hourly data with daily and weekly seasonalities. The final 48 observations are held out to compare the broader automatic search with a cheaper restricted search.\n",
+    "\n",
+    "The timing cells call `fit` and `predict` directly so the fitted model metadata, such as `model_` and `selected_orders_`, is available for inspection.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import time\n",
+    "\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "\n",
+    "from statsforecast import StatsForecast\n",
+    "from statsforecast.models import AutoMSARIMAX\n",
+    "\n",
+    "\n",
+    "def simulate_two_seasonal(n=240, seed=42):\n",
+    "    rng = np.random.default_rng(seed)\n",
+    "    t = np.arange(n)\n",
+    "    daily = 3.0 * np.sin(2 * np.pi * t / 24)\n",
+    "    weekly = 1.8 * np.cos(2 * np.pi * t / 168)\n",
+    "    trend = 0.01 * t\n",
+    "    noise = rng.normal(scale=0.35, size=n)\n",
+    "    return 20 + trend + daily + weekly + noise\n",
+    "\n",
+    "\n",
+    "h = 48\n",
+    "y = simulate_two_seasonal()\n",
+    "df = pd.DataFrame(\n",
+    "    {\n",
+    "        \"unique_id\": \"series_1\",\n",
+    "        \"ds\": pd.date_range(\"2024-01-01\", periods=len(y), freq=\"h\"),\n",
+    "        \"y\": y,\n",
+    "    }\n",
+    ")\n",
+    "train_df = df.iloc[:-h]\n",
+    "train_y = train_df[\"y\"].to_numpy()\n",
+    "test = df[\"y\"].to_numpy()[-h:]\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Fuller Automatic Search\n",
+    "\n",
+    "This configuration keeps the candidate space modest, but it still fits every valid combination inside the specified bounds. With `max_ar_order=[1, 1, 1]`, `max_i_order=[1, 1, 0]`, `max_ma_order=[0, 0, 0]`, and `include_constant=False`, the search evaluates 31 candidate SARIMAX models.\n",
+    "\n",
+    "The number of candidates grows multiplicatively, so adding MA terms or constants can make the same example much slower.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "full_model = AutoMSARIMAX(\n",
+    "    lags=[1, 24, 168],\n",
+    "    max_ar_order=[1, 1, 1],\n",
+    "    max_i_order=[1, 1, 0],\n",
+    "    max_ma_order=[0, 0, 0],\n",
+    "    include_constant=False,\n",
+    "    n_jobs=-1,\n",
+    ")\n",
+    "\n",
+    "start = time.perf_counter()\n",
+    "full_model.fit(train_y)\n",
+    "full_fcst = full_model.predict(h=h, level=[80, 95])\n",
+    "full_elapsed = time.perf_counter() - start\n",
+    "full_rmse = np.sqrt(np.mean((full_fcst[\"mean\"] - test) ** 2))\n",
+    "\n",
+    "print(\"selection:\", full_model.model_[\"selection\"])\n",
+    "print(\"selected orders:\", full_model.selected_orders_)\n",
+    "print(\"elapsed seconds:\", round(full_elapsed, 2))\n",
+    "print(\"48-step RMSE:\", round(full_rmse, 4))\n",
+    "pd.DataFrame(full_fcst).head()\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Cheaper Restricted Search\n",
+    "\n",
+    "A practical first pass is to remove combinations that are unlikely to pay off. The version below still models daily and weekly seasonality, but it fixes non-seasonal AR and differencing to zero, keeps MA terms off, and keeps constants off. This reduces the search from 31 candidate fits to 7.\n",
+    "\n",
+    "Use this kind of restricted search when the broad auto model is too slow, then widen one dimension at a time only if diagnostics or accuracy justify it.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cheap_model = AutoMSARIMAX(\n",
+    "    lags=[1, 24, 168],\n",
+    "    max_ar_order=[0, 1, 1],\n",
+    "    max_i_order=[0, 1, 0],\n",
+    "    max_ma_order=[0, 0, 0],\n",
+    "    include_constant=False,\n",
+    "    n_jobs=-1,\n",
+    ")\n",
+    "\n",
+    "start = time.perf_counter()\n",
+    "cheap_model.fit(train_y)\n",
+    "cheap_fcst = cheap_model.predict(h=h, level=[80, 95])\n",
+    "cheap_elapsed = time.perf_counter() - start\n",
+    "cheap_rmse = np.sqrt(np.mean((cheap_fcst[\"mean\"] - test) ** 2))\n",
+    "\n",
+    "print(\"selection:\", cheap_model.model_[\"selection\"])\n",
+    "print(\"selected orders:\", cheap_model.selected_orders_)\n",
+    "print(\"elapsed seconds:\", round(cheap_elapsed, 2))\n",
+    "print(\"48-step RMSE:\", round(cheap_rmse, 4))\n",
+    "pd.DataFrame(cheap_fcst).head()\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "pygments_lexer": "ipython3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

nbs/docs/models/MSARIMAX.ipynb

@@ -0,0 +1,90 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#| echo: false\n",
+    "import warnings\n",
+    "warnings.filterwarnings(\"ignore\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# MSARIMAX Model\n",
+    "\n",
+    "> Fixed-order multiple-seasonal SARIMAX model for time series with more than one seasonal pattern.\n",
+    "\n",
+    "`MSARIMAX` expands AR, differencing, and MA orders across multiple seasonal lag levels. It supports exogenous regressors, prediction intervals, and in-sample fitted values through the standard StatsForecast model interface."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Example\n",
+    "\n",
+    "The example below simulates hourly data with daily and weekly seasonalities, then fits a fixed multiple-seasonal SARIMAX specification."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "\n",
+    "from statsforecast import StatsForecast\n",
+    "from statsforecast.models import MSARIMAX\n",
+    "\n",
+    "\n",
+    "def simulate_two_seasonal(n=240, seed=42):\n",
+    "    rng = np.random.default_rng(seed)\n",
+    "    t = np.arange(n)\n",
+    "    daily = 3.0 * np.sin(2 * np.pi * t / 24)\n",
+    "    weekly = 1.8 * np.cos(2 * np.pi * t / 168)\n",
+    "    trend = 0.01 * t\n",
+    "    noise = rng.normal(scale=0.35, size=n)\n",
+    "    return 20 + trend + daily + weekly + noise\n",
+    "\n",
+    "\n",
+    "y = simulate_two_seasonal()\n",
+    "df = pd.DataFrame(\n",
+    "    {\n",
+    "        \"unique_id\": \"series_1\",\n",
+    "        \"ds\": pd.date_range(\"2024-01-01\", periods=len(y), freq=\"h\"),\n",
+    "        \"y\": y,\n",
+    "    }\n",
+    ")\n",
+    "\n",
+    "model = MSARIMAX(\n",
+    "    lags=[1, 24, 168],\n",
+    "    ar_order=[0, 0, 1],\n",
+    "    i_order=[0, 1, 0],\n",
+    "    ma_order=[0, 0, 0],\n",
+    ")\n",
+    "sf = StatsForecast(models=[model], freq=\"h\")\n",
+    "sf.forecast(df=df, h=48, level=[80, 95]).head()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "pygments_lexer": "ipython3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}