From 4fefce8b3c55d92ce0b5888135ec1f5db0d4c4e7 Mon Sep 17 00:00:00 2001
From: gdeleva <giacomo.deleva@gmail.com>
Date: Mon, 9 Jun 2025 12:08:53 +0900
Subject: [PATCH 01/16] Add num_steps_per_season parameter in TimeSeasonality

---
 pymc_extras/statespace/models/structural.py | 22 +++++++++++++++++----
 1 file changed, 18 insertions(+), 4 deletions(-)

diff --git a/pymc_extras/statespace/models/structural.py b/pymc_extras/statespace/models/structural.py
index a982366c3..e52434168 100644
--- a/pymc_extras/statespace/models/structural.py
+++ b/pymc_extras/statespace/models/structural.py
@@ -1068,7 +1068,11 @@ class TimeSeasonality(Component):
     ----------
     season_length: int
         The number of periods in a single seasonal cycle, e.g. 12 for monthly data with annual seasonal pattern, 7 for
-        daily data with weekly seasonal pattern, etc.
+        daily data with weekly seasonal pattern, etc. It must be greater than one.
+
+    num_steps_per_season: int, default 1
+        Number of time steps between successive applications of the same seasonal position (state).
+        This determines how long each seasonal effect is held constant before moving to the next.
 
     innovations: bool, default True
         Whether to include stochastic innovations in the strength of the seasonal effect
@@ -1094,14 +1098,24 @@ class TimeSeasonality(Component):
     -----
     A seasonal effect is any pattern that repeats every fixed interval. Although there are many possible ways to
     model seasonal effects, the implementation used here is the one described by [1] as the "canonical" time domain
-    representation. The seasonal component can be expressed:
+    representation. Indexing the seasonal component as
+
+    .. math::
+        \underbrace{\gamma_0, \gamma_0, \ldots, \gamma_0}_{r\ \text{times}},
+        \underbrace{\gamma_1, \gamma_1, \ldots, \gamma_1}_{r\ \text{times}},
+        \ldots,
+        \underbrace{\gamma_{s-1}, \gamma_{s-1}, \ldots, \gamma_{s-1}}_{r\ \text{times}},
+        \ldots,
+
+    where :math:`s` is the ``seasonal_length`` parameter and :math:`r` is the ``num_steps_per_season`` parameter,
+    the seasonal component can be then expressed:
 
     .. math::
         \gamma_t = -\sum_{i=1}^{s-1} \gamma_{t-i} + \omega_t, \quad \omega_t \sim N(0, \sigma_\gamma)
 
-    Where :math:`s` is the ``seasonal_length`` parameter and :math:`\omega_t` is the (optional) stochastic innovation.
+    where :math:`\omega_t` is the (optional) stochastic innovation.
     To give interpretation to the :math:`\gamma` terms, it is helpful to work  through the algebra for a simple
-    example. Let :math:`s=4`, and omit the shock term. Define initial conditions :math:`\gamma_0, \gamma_{-1},
+    example. Let :math:`s=4`, :math:`r=1`, and omit the shock term. Define initial conditions :math:`\gamma_0, \gamma_{-1},
     \gamma_{-2}`. The value of the seasonal component for the first 5 timesteps will be:
 
     .. math::

From 527520cbd93207342defbaaa6285efeda7dacdbb Mon Sep 17 00:00:00 2001
From: gdeleva <giacomo.deleva@gmail.com>
Date: Fri, 27 Jun 2025 19:01:26 +0900
Subject: [PATCH 02/16] Use floor sign to write the seasonal component

---
 pymc_extras/statespace/models/structural.py | 29 +++++++++++----------
 1 file changed, 15 insertions(+), 14 deletions(-)

diff --git a/pymc_extras/statespace/models/structural.py b/pymc_extras/statespace/models/structural.py
index e52434168..cee201cb2 100644
--- a/pymc_extras/statespace/models/structural.py
+++ b/pymc_extras/statespace/models/structural.py
@@ -1070,7 +1070,7 @@ class TimeSeasonality(Component):
         The number of periods in a single seasonal cycle, e.g. 12 for monthly data with annual seasonal pattern, 7 for
         daily data with weekly seasonal pattern, etc. It must be greater than one.
 
-    num_steps_per_season: int, default 1
+    duration: int, default 1
         Number of time steps between successive applications of the same seasonal position (state).
         This determines how long each seasonal effect is held constant before moving to the next.
 
@@ -1098,25 +1098,24 @@ class TimeSeasonality(Component):
     -----
     A seasonal effect is any pattern that repeats every fixed interval. Although there are many possible ways to
     model seasonal effects, the implementation used here is the one described by [1] as the "canonical" time domain
-    representation. Indexing the seasonal component as
+    representation. Given :math:`s` initial states
 
     .. math::
-        \underbrace{\gamma_0, \gamma_0, \ldots, \gamma_0}_{r\ \text{times}},
-        \underbrace{\gamma_1, \gamma_1, \ldots, \gamma_1}_{r\ \text{times}},
-        \ldots,
-        \underbrace{\gamma_{s-1}, \gamma_{s-1}, \ldots, \gamma_{s-1}}_{r\ \text{times}},
-        \ldots,
+        \tilde{\gamma}_{0}, \tilde{\gamma}_{1}, \ldots, \tilde{\gamma}_{s-1},
 
-    where :math:`s` is the ``seasonal_length`` parameter and :math:`r` is the ``num_steps_per_season`` parameter,
-    the seasonal component can be then expressed:
+    where :math:`s` is the ``seasonal_length`` parameter, the full seasonal component can be expressed:
 
     .. math::
-        \gamma_t = -\sum_{i=1}^{s-1} \gamma_{t-i} + \omega_t, \quad \omega_t \sim N(0, \sigma_\gamma)
+        \begin{align}
+            \gamma_t &= \tilde{\gamma}_{k_t}, \quad \text{where} \quad k_t = \left\lfloor \frac{t}{d} \right\rfloor \bmod s \\
+            \tilde{\gamma}_k &= -\sum_{i=1}^{s-1} \tilde{\gamma}_{k - i} + \omega_k, \quad \omega_k \sim \mathcal{N}(0, \sigma)
+        \end{align}
+
+    where :math:`d` is the ``duration`` parameter and :math:`\omega_t` is the (optional) stochastic innovation.
 
-    where :math:`\omega_t` is the (optional) stochastic innovation.
     To give interpretation to the :math:`\gamma` terms, it is helpful to work  through the algebra for a simple
-    example. Let :math:`s=4`, :math:`r=1`, and omit the shock term. Define initial conditions :math:`\gamma_0, \gamma_{-1},
-    \gamma_{-2}`. The value of the seasonal component for the first 5 timesteps will be:
+    example. Let :math:`s=4`, :math:`d=1`, and omit the shock term. Define initial conditions :math:`\tilde{\gamma}_0, \tilde{\gamma}_{1},
+    \tilde{\gamma}_{2}`. The value of the seasonal component for the first 5 timesteps will be:
 
     .. math::
         \begin{align}
@@ -1193,13 +1192,14 @@ class TimeSeasonality(Component):
     def __init__(
         self,
         season_length: int,
+        duration: int = 1,
         innovations: bool = True,
         name: str | None = None,
         state_names: list | None = None,
         remove_first_state: bool = True,
     ):
         if name is None:
-            name = f"Seasonal[s={season_length}]"
+            name = f"Seasonal[s={season_length}, d={duration}]"
         if state_names is None:
             state_names = [f"{name}_{i}" for i in range(season_length)]
         else:
@@ -1209,6 +1209,7 @@ def __init__(
                 )
             state_names = state_names.copy()
         self.innovations = innovations
+        self.duration = duration
         self.remove_first_state = remove_first_state
 
         if self.remove_first_state:

From 2b7eede174819099c87a32d94da30230dbd4574d Mon Sep 17 00:00:00 2001
From: gdeleva <giacomo.deleva@gmail.com>
Date: Mon, 14 Jul 2025 18:20:16 +0900
Subject: [PATCH 03/16] Models Definition

---
 pymc_extras/statespace/models/structural.py | 84 +++++++++++++++++----
 1 file changed, 68 insertions(+), 16 deletions(-)

diff --git a/pymc_extras/statespace/models/structural.py b/pymc_extras/statespace/models/structural.py
index cee201cb2..8f1debf2a 100644
--- a/pymc_extras/statespace/models/structural.py
+++ b/pymc_extras/statespace/models/structural.py
@@ -1071,48 +1071,98 @@ class TimeSeasonality(Component):
         daily data with weekly seasonal pattern, etc. It must be greater than one.
 
     duration: int, default 1
-        Number of time steps between successive applications of the same seasonal position (state).
-        This determines how long each seasonal effect is held constant before moving to the next.
+        Number of time steps for each seasonal period.
+        This determines how long each seasonal period is held constant before moving to the next.
 
     innovations: bool, default True
         Whether to include stochastic innovations in the strength of the seasonal effect
 
     name: str, default None
         A name for this seasonal component. Used to label dimensions and coordinates. Useful when multiple seasonal
-        components are included in the same model. Default is ``f"Seasonal[s={season_length}]"``
+        components are included in the same model. Default is ``f"Seasonal[s={season_length}, d={duration}]"``
 
     state_names: list of str, default None
         List of strings for seasonal effect labels. If provided, it must be of length ``season_length``. An example
         would be ``state_names = ['Mon', 'Tue', 'Wed', 'Thur', 'Fri', 'Sat', 'Sun']`` when data is daily with a weekly
         seasonal pattern (``season_length = 7``).
 
-        If None, states will be numbered ``[State_0, ..., State_s]``
+        If None, states will be numbered ``[State_0, ..., State_s-1]``
 
     remove_first_state: bool, default True
-        If True, the first state will be removed from the model. This is done because there are only n-1 degrees of
+        If True, the first state will be removed from the model. This is done because there are only ``season_length-1`` degrees of
         freedom in the seasonal component, and one state is not identified. If False, the first state will be
         included in the model, but it will not be identified -- you will need to handle this in the priors (e.g. with
         ZeroSumNormal).
 
     Notes
     -----
-    A seasonal effect is any pattern that repeats every fixed interval. Although there are many possible ways to
-    model seasonal effects, the implementation used here is the one described by [1] as the "canonical" time domain
-    representation. Given :math:`s` initial states
+    A seasonal effect is any pattern that repeats at fixed intervals. There are several ways to model such effects;
+    here, we present two models that are straightforward extensions of those described in [1].
+
+    First model (``remove_first_state=True``)
+
+    In this model, the state vector is defined as:
 
     .. math::
-        \tilde{\gamma}_{0}, \tilde{\gamma}_{1}, \ldots, \tilde{\gamma}_{s-1},
+        \alpha_t :=(\gamma_t, \ldots, \gamma_{t-d(s-1)+1}), \quad t \ge 1.
+
+    This vector has length :math:`d(s-1)`, where:
 
-    where :math:`s` is the ``seasonal_length`` parameter, the full seasonal component can be expressed:
+    - :math:`s` is the ``seasonal_length`` parameter, and
+    - :math:`d` is the ``duration`` parameter.
+
+    The components of the initial vector :math:`\alpha_{1}` are given by
 
     .. math::
-        \begin{align}
-            \gamma_t &= \tilde{\gamma}_{k_t}, \quad \text{where} \quad k_t = \left\lfloor \frac{t}{d} \right\rfloor \bmod s \\
-            \tilde{\gamma}_k &= -\sum_{i=1}^{s-1} \tilde{\gamma}_{k - i} + \omega_k, \quad \omega_k \sim \mathcal{N}(0, \sigma)
-        \end{align}
+        \gamma_{1-l} := \tilde{\gamma}_{1+k_l}, \quad \text{where} \quad k_l := \left\lfloor \frac{l}{d} \right\rfloor \bmod s \quad \text{and} \quad l=0,\ldots, d(s-1)-1.
+
+    Here, the values
+
+    .. math::
+        \tilde{\gamma}_{1}, \ldots, \tilde{\gamma}_{s-1},
+
+    represent the initial seasonal states. The transition matrix of this model is the :math:`d(s-1) \times d(s-1)` matrix
+
+    .. math::
+        \begin{bmatrix}
+            -\mathbf{1}_d & -\mathbf{1}_d & \cdots & -\mathbf{1}_d \\
+            \mathbf{0}_d & \mathbf{1}_d & \cdots & \mathbf{0}_d \\
+            \vdots & \vdots & \ddots & \vdots \\
+            \mathbf{0}_d & \mathbf{0}_d & \cdots & \mathbf{1}_d
+        \end{bmatrix}
+
+    where :math:`\mathbf{1}_d` and  :math:`\mathbf{0}_d` denote the :math:`d \times d` identity and null matrices, respectively.
 
-    where :math:`d` is the ``duration`` parameter and :math:`\omega_t` is the (optional) stochastic innovation.
+    Second model (``remove_first_state=False``)
 
+    In contrast, the state vector in the second model is defined as:
+
+    .. math::
+        \alpha_t=(\gamma_t, \ldots, \gamma_{t-ds+1}), \quad t \ge 1.
+
+    This vector has length :math:`ds`. The components of the initial state vector :math:`\alpha_{1}` are defined similarly:
+
+    .. math::
+        \gamma_{1-l} := \tilde{\gamma}_{1+k_l}, \quad \text{where} \quad k_l := \left\lfloor \frac{l}{d} \right\rfloor \bmod s \quad \text{and} \quad l=0,\ldots, ds-1.
+
+    In this case, the initial seasonal states are required to satisfy the following condition:
+
+    .. math::
+        \sum_{i=1}^{s} \tilde{\gamma}_{i} = 0.
+
+    The transition matrix of this model is the following :math:`ds \times ds` circulant matrix:
+
+    .. math::
+        \begin{bmatrix}
+            0 & 1 & 0 & \cdots & 0 \\
+            0 & 0 & 1 & \cdots & 0 \\
+            \vdots & \vdots & \ddots & \ddots & \vdots \\
+            0 & 0 & \cdots & 0 & 1 \\
+            1 & 0 & \cdots & 0 & 0
+        \end{bmatrix}
+
+    Examples
+    --------
     To give interpretation to the :math:`\gamma` terms, it is helpful to work  through the algebra for a simple
     example. Let :math:`s=4`, :math:`d=1`, and omit the shock term. Define initial conditions :math:`\tilde{\gamma}_0, \tilde{\gamma}_{1},
     \tilde{\gamma}_{2}`. The value of the seasonal component for the first 5 timesteps will be:
@@ -1198,6 +1248,8 @@ def __init__(
         state_names: list | None = None,
         remove_first_state: bool = True,
     ):
+        if season_length <= 1:
+            raise ValueError(f"season_length must be greater than 1, got {season_length}")
         if name is None:
             name = f"Seasonal[s={season_length}, d={duration}]"
         if state_names is None:
@@ -1218,7 +1270,7 @@ def __init__(
             # TODO: Can this be stashed and reconstructed automatically somehow?
             state_names.pop(0)
 
-        k_states = season_length - int(self.remove_first_state)
+        k_states = season_length * duration - int(self.remove_first_state) * duration
 
         super().__init__(
             name=name,

From d7ff1d4917531a0354e6449b79116565572717d3 Mon Sep 17 00:00:00 2001
From: gdeleva <giacomo.deleva@gmail.com>
Date: Tue, 15 Jul 2025 22:50:13 +0900
Subject: [PATCH 04/16] Edit symbolic graph

---
 pymc_extras/statespace/models/structural.py | 47 +++++++++++++++------
 1 file changed, 33 insertions(+), 14 deletions(-)

diff --git a/pymc_extras/statespace/models/structural.py b/pymc_extras/statespace/models/structural.py
index 8f1debf2a..6ccdcb436 100644
--- a/pymc_extras/statespace/models/structural.py
+++ b/pymc_extras/statespace/models/structural.py
@@ -1099,7 +1099,7 @@ class TimeSeasonality(Component):
     A seasonal effect is any pattern that repeats at fixed intervals. There are several ways to model such effects;
     here, we present two models that are straightforward extensions of those described in [1].
 
-    First model (``remove_first_state=True``)
+    **First model** (``remove_first_state=True``)
 
     In this model, the state vector is defined as:
 
@@ -1125,15 +1125,16 @@ class TimeSeasonality(Component):
 
     .. math::
         \begin{bmatrix}
-            -\mathbf{1}_d & -\mathbf{1}_d & \cdots & -\mathbf{1}_d \\
-            \mathbf{0}_d & \mathbf{1}_d & \cdots & \mathbf{0}_d \\
+            -\mathbf{1}_d & -\mathbf{1}_d & \cdots & -\mathbf{1}_d & -\mathbf{1}_d \\
+            \mathbf{1}_d & \mathbf{0}_d & \cdots & \mathbf{0}_d & \mathbf{0}_d \\
+            \mathbf{0}_d & \mathbf{1}_d & \cdots & \mathbf{0}_d & \mathbf{0}_d \\
             \vdots & \vdots & \ddots & \vdots \\
-            \mathbf{0}_d & \mathbf{0}_d & \cdots & \mathbf{1}_d
+            \mathbf{0}_d & \mathbf{0}_d & \cdots & \mathbf{1}_d & \mathbf{0}_d
         \end{bmatrix}
 
     where :math:`\mathbf{1}_d` and  :math:`\mathbf{0}_d` denote the :math:`d \times d` identity and null matrices, respectively.
 
-    Second model (``remove_first_state=False``)
+    **Second model** (``remove_first_state=False``)
 
     In contrast, the state vector in the second model is defined as:
 
@@ -1161,9 +1162,7 @@ class TimeSeasonality(Component):
             1 & 0 & \cdots & 0 & 0
         \end{bmatrix}
 
-    Examples
-    --------
-    To give interpretation to the :math:`\gamma` terms, it is helpful to work  through the algebra for a simple
+    To give interpretation to the :math:`\gamma` terms, it is helpful to work through the algebra for a simple
     example. Let :math:`s=4`, :math:`d=1`, and omit the shock term. Define initial conditions :math:`\tilde{\gamma}_0, \tilde{\gamma}_{1},
     \tilde{\gamma}_{2}`. The value of the seasonal component for the first 5 timesteps will be:
 
@@ -1253,11 +1252,11 @@ def __init__(
         if name is None:
             name = f"Seasonal[s={season_length}, d={duration}]"
         if state_names is None:
-            state_names = [f"{name}_{i}" for i in range(season_length)]
+            state_names = [f"{name}_{i}_{j}" for i in range(season_length) for j in range(duration)]
         else:
-            if len(state_names) != season_length:
+            if len(state_names) != season_length * duration:
                 raise ValueError(
-                    f"state_names must be a list of length season_length, got {len(state_names)}"
+                    f"state_names must be a list of length season_length*duration, got {len(state_names)}"
                 )
             state_names = state_names.copy()
         self.innovations = innovations
@@ -1268,7 +1267,7 @@ def __init__(
             # In traditional models, the first state isn't identified, so we can help out the user by automatically
             # discarding it.
             # TODO: Can this be stashed and reconstructed automatically somehow?
-            state_names.pop(0)
+            state_names = state_names[duration:]
 
         k_states = season_length * duration - int(self.remove_first_state) * duration
 
@@ -1308,8 +1307,28 @@ def make_symbolic_graph(self) -> None:
         if self.remove_first_state:
             # In this case, parameters are normalized to sum to zero, so the current state is the negative sum of
             # all previous states.
-            T = np.eye(self.k_states, k=-1)
-            T[0, :] = -1
+            one_d = np.ones((self.duration, self.duration))
+            zero_d = np.zeros((self.duration, self.duration))
+            id_d = np.eye(self.duration)
+
+            blocks = []
+
+            # First row: all -1_d blocks
+            first_row = [-one_d for _ in range(self.season_length - 1)]
+            blocks.append(first_row)
+
+            # Rows 2 to season_length-1: shifted identity blocks
+            for i in range(self.season_length - 2):
+                row = []
+                for j in range(self.season_length - 1):
+                    if j == i:
+                        row.append(id_d)
+                    else:
+                        row.append(zero_d)
+                blocks.append(row)
+
+            # Stack blocks
+            T = np.block(blocks)
         else:
             # In this case we assume the user to be responsible for ensuring the states sum to zero, so T is just a
             # circulant matrix that cycles between the states.

From 639aa14c5b65c827fba7ef2a93c041140aa150d1 Mon Sep 17 00:00:00 2001
From: gdeleva <giacomo.deleva@gmail.com>
Date: Wed, 16 Jul 2025 10:03:34 +0900
Subject: [PATCH 05/16] Correct transition matrix

---
 pymc_extras/statespace/models/structural.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/pymc_extras/statespace/models/structural.py b/pymc_extras/statespace/models/structural.py
index 6ccdcb436..a02e55da4 100644
--- a/pymc_extras/statespace/models/structural.py
+++ b/pymc_extras/statespace/models/structural.py
@@ -1307,14 +1307,13 @@ def make_symbolic_graph(self) -> None:
         if self.remove_first_state:
             # In this case, parameters are normalized to sum to zero, so the current state is the negative sum of
             # all previous states.
-            one_d = np.ones((self.duration, self.duration))
             zero_d = np.zeros((self.duration, self.duration))
             id_d = np.eye(self.duration)
 
             blocks = []
 
             # First row: all -1_d blocks
-            first_row = [-one_d for _ in range(self.season_length - 1)]
+            first_row = [-id_d for _ in range(self.season_length - 1)]
             blocks.append(first_row)
 
             # Rows 2 to season_length-1: shifted identity blocks

From c1937b4b59af877f17f2a1ec3069e10225e44804 Mon Sep 17 00:00:00 2001
From: gdeleva <giacomo.deleva@gmail.com>
Date: Wed, 16 Jul 2025 12:12:55 +0900
Subject: [PATCH 06/16] Correct documentation

---
 pymc_extras/statespace/models/structural.py | 23 +++++++++++----------
 1 file changed, 12 insertions(+), 11 deletions(-)

diff --git a/pymc_extras/statespace/models/structural.py b/pymc_extras/statespace/models/structural.py
index a02e55da4..69ffc9edc 100644
--- a/pymc_extras/statespace/models/structural.py
+++ b/pymc_extras/statespace/models/structural.py
@@ -1104,22 +1104,22 @@ class TimeSeasonality(Component):
     In this model, the state vector is defined as:
 
     .. math::
-        \alpha_t :=(\gamma_t, \ldots, \gamma_{t-d(s-1)+1}), \quad t \ge 1.
+        \alpha_t :=(\gamma_t, \ldots, \gamma_{t-d(s-1)+1}), \quad t \ge 0.
 
     This vector has length :math:`d(s-1)`, where:
 
     - :math:`s` is the ``seasonal_length`` parameter, and
     - :math:`d` is the ``duration`` parameter.
 
-    The components of the initial vector :math:`\alpha_{1}` are given by
+    The components of the initial vector :math:`\alpha_{0}` are given by
 
     .. math::
-        \gamma_{1-l} := \tilde{\gamma}_{1+k_l}, \quad \text{where} \quad k_l := \left\lfloor \frac{l}{d} \right\rfloor \bmod s \quad \text{and} \quad l=0,\ldots, d(s-1)-1.
+        \gamma_{-l} := \tilde{\gamma}_{k_l}, \quad \text{where} \quad k_l := \left\lfloor \frac{l}{d} \right\rfloor \bmod s \quad \text{and} \quad l=0,\ldots, d(s-1)-1.
 
     Here, the values
 
     .. math::
-        \tilde{\gamma}_{1}, \ldots, \tilde{\gamma}_{s-1},
+        \tilde{\gamma}_{0}, \ldots, \tilde{\gamma}_{s-2},
 
     represent the initial seasonal states. The transition matrix of this model is the :math:`d(s-1) \times d(s-1)` matrix
 
@@ -1139,17 +1139,17 @@ class TimeSeasonality(Component):
     In contrast, the state vector in the second model is defined as:
 
     .. math::
-        \alpha_t=(\gamma_t, \ldots, \gamma_{t-ds+1}), \quad t \ge 1.
+        \alpha_t=(\gamma_t, \ldots, \gamma_{t-ds+1}), \quad t \ge 0.
 
-    This vector has length :math:`ds`. The components of the initial state vector :math:`\alpha_{1}` are defined similarly:
+    This vector has length :math:`ds`. The components of the initial state vector :math:`\alpha_{0}` are defined similarly:
 
     .. math::
-        \gamma_{1-l} := \tilde{\gamma}_{1+k_l}, \quad \text{where} \quad k_l := \left\lfloor \frac{l}{d} \right\rfloor \bmod s \quad \text{and} \quad l=0,\ldots, ds-1.
+        \gamma_{-l} := \tilde{\gamma}_{k_l}, \quad \text{where} \quad k_l := \left\lfloor \frac{l}{d} \right\rfloor \bmod s \quad \text{and} \quad l=0,\ldots, ds-1.
 
-    In this case, the initial seasonal states are required to satisfy the following condition:
+    In this case, the initial seasonal states :math:`\tilde{\gamma}_{0}, \ldots, \tilde{\gamma}_{s-1}` are required to satisfy the following condition:
 
     .. math::
-        \sum_{i=1}^{s} \tilde{\gamma}_{i} = 0.
+        \sum_{i=0}^{s-1} \tilde{\gamma}_{i} = 0.
 
     The transition matrix of this model is the following :math:`ds \times ds` circulant matrix:
 
@@ -1163,8 +1163,9 @@ class TimeSeasonality(Component):
         \end{bmatrix}
 
     To give interpretation to the :math:`\gamma` terms, it is helpful to work through the algebra for a simple
-    example. Let :math:`s=4`, :math:`d=1`, and omit the shock term. Define initial conditions :math:`\tilde{\gamma}_0, \tilde{\gamma}_{1},
-    \tilde{\gamma}_{2}`. The value of the seasonal component for the first 5 timesteps will be:
+    example. Let :math:`s=4`, :math:`d=1`, ``remove_first_state=True``, and omit the shock term. Then, we have
+    :math:`\gamma_{-i} = \tilde{\gamma}_{-i}`, for :math:`i=-2,\ldots, 0` and the value of the seasonal component
+    for the first 5 timesteps will be:
 
     .. math::
         \begin{align}

From 13ffcddafc741a9f5d680ea90a91900bbd6a8863 Mon Sep 17 00:00:00 2001
From: gdeleva <giacomo.deleva@gmail.com>
Date: Wed, 16 Jul 2025 14:40:28 +0900
Subject: [PATCH 07/16] Add second examples in the documentation

---
 pymc_extras/statespace/models/structural.py | 34 +++++++++++++++++++--
 1 file changed, 31 insertions(+), 3 deletions(-)

diff --git a/pymc_extras/statespace/models/structural.py b/pymc_extras/statespace/models/structural.py
index 69ffc9edc..63b260279 100644
--- a/pymc_extras/statespace/models/structural.py
+++ b/pymc_extras/statespace/models/structural.py
@@ -1199,10 +1199,38 @@ class TimeSeasonality(Component):
     And so on. So for interpretation, the ``season_length - 1`` initial states are, when reversed, the coefficients
     associated with ``state_names[1:]``.
 
+    In the next example, we set :math:`s=3`, :math:`d=2`, ``remove_first_state=True``, and omit the shock term.
+    By definition, the initial vector :math:`\alpha_{0}` is
+
+    .. math::
+        \alpha_0=(\tilde{\gamma}_{0}, \tilde{\gamma}_{0}, \tilde{\gamma}_{-1}, \tilde{\gamma}_{-1})
+
+    and the transition matrix is
+
+    .. math::
+        \begin{bmatrix}
+            -1 &  0 & -1 &  0 \\
+             0 & -1 &  0 & -1 \\
+             1 &  0 &  0 &  0 \\
+             0 &  1 &  0 &  0 \\
+        \end{bmatrix}
+
+    It is easy to verify that:
+
+    .. math::
+        \begin{align}
+            \gamma_1 &= -\tilde{\gamma}_0 - \tilde{\gamma}_{-1}\\
+            \gamma_2 &= -(-\tilde{\gamma}_0 - \tilde{\gamma}_{-1})-\tilde{\gamma}_0\\
+                     &= \tilde{\gamma}_{-1}\\
+            \gamma_3 &= -\tilde{\gamma}_{-1} +(\tilde{\gamma}_0 + \tilde{\gamma}_{-1})\\
+                     &= \tilde{\gamma}_{0}\\
+            \gamma_4 &= -\tilde{\gamma}_0 - \tilde{\gamma}_{-1}.\\
+        \end{align}
+
     .. warning::
-        Although the ``state_names`` argument expects a list of length ``season_length``, only ``state_names[1:]``
-        will be saved as model dimensions, since the 1st coefficient is not identified (it is defined as
-        :math:`-\sum_{i=1}^{s} \gamma_{t-i}`).
+        Although the ``state_names`` argument expects a list of length ``season_length`` times ``duration``,
+        only ``state_names[1:]`` will be saved as model dimensions, since the first coefficient is not identified
+        (it is defined as :math:`-\sum_{i=1}^{s-1} \tilde{\gamma}_{-i}`).
 
     Examples
     --------

From 0868f90ee8a7cc9360aa7aea7e509a31e6786019 Mon Sep 17 00:00:00 2001
From: gdeleva <giacomo.deleva@gmail.com>
Date: Thu, 17 Jul 2025 10:25:24 +0900
Subject: [PATCH 08/16] Fix bugs

---
 pymc_extras/statespace/models/structural.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/pymc_extras/statespace/models/structural.py b/pymc_extras/statespace/models/structural.py
index 63b260279..9c8b5a108 100644
--- a/pymc_extras/statespace/models/structural.py
+++ b/pymc_extras/statespace/models/structural.py
@@ -1281,7 +1281,12 @@ def __init__(
         if name is None:
             name = f"Seasonal[s={season_length}, d={duration}]"
         if state_names is None:
-            state_names = [f"{name}_{i}_{j}" for i in range(season_length) for j in range(duration)]
+            if duration > 1:
+                state_names = [
+                    f"{name}_{i}_{j}" for i in range(season_length) for j in range(duration)
+                ]
+            else:
+                state_names = [f"{name}_{i}" for i in range(season_length)]
         else:
             if len(state_names) != season_length * duration:
                 raise ValueError(
@@ -1291,6 +1296,7 @@ def __init__(
         self.innovations = innovations
         self.duration = duration
         self.remove_first_state = remove_first_state
+        self.season_length = season_length
 
         if self.remove_first_state:
             # In traditional models, the first state isn't identified, so we can help out the user by automatically

From 07d94b61c52bd6799fd197aeaf4cd3165fd6f951 Mon Sep 17 00:00:00 2001
From: gdeleva <giacomo.deleva@gmail.com>
Date: Thu, 17 Jul 2025 11:04:46 +0900
Subject: [PATCH 09/16] Merge with main

---
 .../structural/components/seasonality.py      | 166 +++++++++++++++---
 1 file changed, 143 insertions(+), 23 deletions(-)

diff --git a/pymc_extras/statespace/models/structural/components/seasonality.py b/pymc_extras/statespace/models/structural/components/seasonality.py
index 1d103fbf0..9df859adc 100644
--- a/pymc_extras/statespace/models/structural/components/seasonality.py
+++ b/pymc_extras/statespace/models/structural/components/seasonality.py
@@ -14,24 +14,28 @@ class TimeSeasonality(Component):
     ----------
     season_length: int
         The number of periods in a single seasonal cycle, e.g. 12 for monthly data with annual seasonal pattern, 7 for
-        daily data with weekly seasonal pattern, etc.
+        daily data with weekly seasonal pattern, etc. It must be greater than one.
+
+    duration: int, default 1
+        Number of time steps for each seasonal period.
+        This determines how long each seasonal period is held constant before moving to the next.
 
     innovations: bool, default True
         Whether to include stochastic innovations in the strength of the seasonal effect
 
     name: str, default None
         A name for this seasonal component. Used to label dimensions and coordinates. Useful when multiple seasonal
-        components are included in the same model. Default is ``f"Seasonal[s={season_length}]"``
+        components are included in the same model. Default is ``f"Seasonal[s={season_length}, d={duration}]"``
 
     state_names: list of str, default None
         List of strings for seasonal effect labels. If provided, it must be of length ``season_length``. An example
         would be ``state_names = ['Mon', 'Tue', 'Wed', 'Thur', 'Fri', 'Sat', 'Sun']`` when data is daily with a weekly
         seasonal pattern (``season_length = 7``).
 
-        If None, states will be numbered ``[State_0, ..., State_s]``
+        If None, states will be numbered ``[State_0, ..., State_s-1]``
 
     remove_first_state: bool, default True
-        If True, the first state will be removed from the model. This is done because there are only n-1 degrees of
+        If True, the first state will be removed from the model. This is done because there are only ``season_length-1`` degrees of
         freedom in the seasonal component, and one state is not identified. If False, the first state will be
         included in the model, but it will not be identified -- you will need to handle this in the priors (e.g. with
         ZeroSumNormal).
@@ -41,17 +45,76 @@ class TimeSeasonality(Component):
 
     Notes
     -----
-    A seasonal effect is any pattern that repeats every fixed interval. Although there are many possible ways to
-    model seasonal effects, the implementation used here is the one described by [1] as the "canonical" time domain
-    representation. The seasonal component can be expressed:
+    A seasonal effect is any pattern that repeats at fixed intervals. There are several ways to model such effects;
+    here, we present two models that are straightforward extensions of those described in [1].
+
+    **First model** (``remove_first_state=True``)
+
+    In this model, the state vector is defined as:
+
+    .. math::
+        \alpha_t :=(\gamma_t, \ldots, \gamma_{t-d(s-1)+1}), \quad t \ge 0.
+
+    This vector has length :math:`d(s-1)`, where:
+
+    - :math:`s` is the ``seasonal_length`` parameter, and
+    - :math:`d` is the ``duration`` parameter.
+
+    The components of the initial vector :math:`\alpha_{0}` are given by
+
+    .. math::
+        \gamma_{-l} := \tilde{\gamma}_{k_l}, \quad \text{where} \quad k_l := \left\lfloor \frac{l}{d} \right\rfloor \bmod s \quad \text{and} \quad l=0,\ldots, d(s-1)-1.
+
+    Here, the values
+
+    .. math::
+        \tilde{\gamma}_{0}, \ldots, \tilde{\gamma}_{s-2},
+
+    represent the initial seasonal states. The transition matrix of this model is the :math:`d(s-1) \times d(s-1)` matrix
+
+    .. math::
+        \begin{bmatrix}
+            -\mathbf{1}_d & -\mathbf{1}_d & \cdots & -\mathbf{1}_d & -\mathbf{1}_d \\
+            \mathbf{1}_d & \mathbf{0}_d & \cdots & \mathbf{0}_d & \mathbf{0}_d \\
+            \mathbf{0}_d & \mathbf{1}_d & \cdots & \mathbf{0}_d & \mathbf{0}_d \\
+            \vdots & \vdots & \ddots & \vdots \\
+            \mathbf{0}_d & \mathbf{0}_d & \cdots & \mathbf{1}_d & \mathbf{0}_d
+        \end{bmatrix}
+
+    where :math:`\mathbf{1}_d` and  :math:`\mathbf{0}_d` denote the :math:`d \times d` identity and null matrices, respectively.
+
+    **Second model** (``remove_first_state=False``)
+
+    In contrast, the state vector in the second model is defined as:
 
     .. math::
-        \gamma_t = -\sum_{i=1}^{s-1} \gamma_{t-i} + \omega_t, \quad \omega_t \sim N(0, \sigma_\gamma)
+        \alpha_t=(\gamma_t, \ldots, \gamma_{t-ds+1}), \quad t \ge 0.
 
-    Where :math:`s` is the ``seasonal_length`` parameter and :math:`\omega_t` is the (optional) stochastic innovation.
-    To give interpretation to the :math:`\gamma` terms, it is helpful to work  through the algebra for a simple
-    example. Let :math:`s=4`, and omit the shock term. Define initial conditions :math:`\gamma_0, \gamma_{-1},
-    \gamma_{-2}`. The value of the seasonal component for the first 5 timesteps will be:
+    This vector has length :math:`ds`. The components of the initial state vector :math:`\alpha_{0}` are defined similarly:
+
+    .. math::
+        \gamma_{-l} := \tilde{\gamma}_{k_l}, \quad \text{where} \quad k_l := \left\lfloor \frac{l}{d} \right\rfloor \bmod s \quad \text{and} \quad l=0,\ldots, ds-1.
+
+    In this case, the initial seasonal states :math:`\tilde{\gamma}_{0}, \ldots, \tilde{\gamma}_{s-1}` are required to satisfy the following condition:
+
+    .. math::
+        \sum_{i=0}^{s-1} \tilde{\gamma}_{i} = 0.
+
+    The transition matrix of this model is the following :math:`ds \times ds` circulant matrix:
+
+    .. math::
+        \begin{bmatrix}
+            0 & 1 & 0 & \cdots & 0 \\
+            0 & 0 & 1 & \cdots & 0 \\
+            \vdots & \vdots & \ddots & \ddots & \vdots \\
+            0 & 0 & \cdots & 0 & 1 \\
+            1 & 0 & \cdots & 0 & 0
+        \end{bmatrix}
+
+    To give interpretation to the :math:`\gamma` terms, it is helpful to work through the algebra for a simple
+    example. Let :math:`s=4`, :math:`d=1`, ``remove_first_state=True``, and omit the shock term. Then, we have
+    :math:`\gamma_{-i} = \tilde{\gamma}_{-i}`, for :math:`i=-2,\ldots, 0` and the value of the seasonal component
+    for the first 5 timesteps will be:
 
     .. math::
         \begin{align}
@@ -85,10 +148,38 @@ class TimeSeasonality(Component):
     And so on. So for interpretation, the ``season_length - 1`` initial states are, when reversed, the coefficients
     associated with ``state_names[1:]``.
 
+    In the next example, we set :math:`s=3`, :math:`d=2`, ``remove_first_state=True``, and omit the shock term.
+    By definition, the initial vector :math:`\alpha_{0}` is
+
+    .. math::
+        \alpha_0=(\tilde{\gamma}_{0}, \tilde{\gamma}_{0}, \tilde{\gamma}_{-1}, \tilde{\gamma}_{-1})
+
+    and the transition matrix is
+
+    .. math::
+        \begin{bmatrix}
+            -1 &  0 & -1 &  0 \\
+             0 & -1 &  0 & -1 \\
+             1 &  0 &  0 &  0 \\
+             0 &  1 &  0 &  0 \\
+        \end{bmatrix}
+
+    It is easy to verify that:
+
+    .. math::
+        \begin{align}
+            \gamma_1 &= -\tilde{\gamma}_0 - \tilde{\gamma}_{-1}\\
+            \gamma_2 &= -(-\tilde{\gamma}_0 - \tilde{\gamma}_{-1})-\tilde{\gamma}_0\\
+                     &= \tilde{\gamma}_{-1}\\
+            \gamma_3 &= -\tilde{\gamma}_{-1} +(\tilde{\gamma}_0 + \tilde{\gamma}_{-1})\\
+                     &= \tilde{\gamma}_{0}\\
+            \gamma_4 &= -\tilde{\gamma}_0 - \tilde{\gamma}_{-1}.\\
+        \end{align}
+
     .. warning::
-        Although the ``state_names`` argument expects a list of length ``season_length``, only ``state_names[1:]``
-        will be saved as model dimensions, since the 1st coefficient is not identified (it is defined as
-        :math:`-\sum_{i=1}^{s} \gamma_{t-i}`).
+        Although the ``state_names`` argument expects a list of length ``season_length`` times ``duration``,
+        only ``state_names[1:]`` will be saved as model dimensions, since the first coefficient is not identified
+        (it is defined as :math:`-\sum_{i=1}^{s-1} \tilde{\gamma}_{-i}`).
 
     Examples
     --------
@@ -137,6 +228,7 @@ class TimeSeasonality(Component):
     def __init__(
         self,
         season_length: int,
+        duration: int = 1,
         innovations: bool = True,
         name: str | None = None,
         state_names: list | None = None,
@@ -146,29 +238,38 @@ def __init__(
         if observed_state_names is None:
             observed_state_names = ["data"]
 
+        if season_length <= 1:
+            raise ValueError(f"season_length must be greater than 1, got {season_length}")
         if name is None:
-            name = f"Seasonal[s={season_length}]"
+            name = f"Seasonal[s={season_length}, d={duration}]"
         if state_names is None:
-            state_names = [f"{name}_{i}" for i in range(season_length)]
+            if duration > 1:
+                state_names = [
+                    f"{name}_{i}_{j}" for i in range(season_length) for j in range(duration)
+                ]
+            else:
+                state_names = [f"{name}_{i}" for i in range(season_length)]
         else:
-            if len(state_names) != season_length:
+            if len(state_names) != season_length * duration:
                 raise ValueError(
-                    f"state_names must be a list of length season_length, got {len(state_names)}"
+                    f"state_names must be a list of length season_length*duration, got {len(state_names)}"
                 )
             state_names = state_names.copy()
 
         self.innovations = innovations
+        self.duration = duration
         self.remove_first_state = remove_first_state
+        self.season_length = season_length
 
         if self.remove_first_state:
             # In traditional models, the first state isn't identified, so we can help out the user by automatically
             # discarding it.
             # TODO: Can this be stashed and reconstructed automatically somehow?
-            state_names.pop(0)
+            state_names = state_names[duration:]
 
         self.provided_state_names = state_names
 
-        k_states = season_length - int(self.remove_first_state)
+        k_states = season_length * duration - int(self.remove_first_state) * duration
         k_endog = len(observed_state_names)
         k_posdef = int(innovations)
 
@@ -236,8 +337,27 @@ def make_symbolic_graph(self) -> None:
         if self.remove_first_state:
             # In this case, parameters are normalized to sum to zero, so the current state is the negative sum of
             # all previous states.
-            T = np.eye(k_states, k=-1)
-            T[0, :] = -1
+            zero_d = np.zeros((self.duration, self.duration))
+            id_d = np.eye(self.duration)
+
+            blocks = []
+
+            # First row: all -1_d blocks
+            first_row = [-id_d for _ in range(self.season_length - 1)]
+            blocks.append(first_row)
+
+            # Rows 2 to season_length-1: shifted identity blocks
+            for i in range(self.season_length - 2):
+                row = []
+                for j in range(self.season_length - 1):
+                    if j == i:
+                        row.append(id_d)
+                    else:
+                        row.append(zero_d)
+                blocks.append(row)
+
+            # Stack blocks
+            T = np.block(blocks)
         else:
             # In this case we assume the user to be responsible for ensuring the states sum to zero, so T is just a
             # circulant matrix that cycles between the states.

From 452f1cb3f63059c908e4bbe1b361ece7b4bdc308 Mon Sep 17 00:00:00 2001
From: gdeleva <giacomo.deleva@gmail.com>
Date: Fri, 18 Jul 2025 10:58:02 +0900
Subject: [PATCH 10/16] Minor fixes

---
 .../models/structural/components/seasonality.py     | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/pymc_extras/statespace/models/structural/components/seasonality.py b/pymc_extras/statespace/models/structural/components/seasonality.py
index 9df859adc..9769ff95a 100644
--- a/pymc_extras/statespace/models/structural/components/seasonality.py
+++ b/pymc_extras/statespace/models/structural/components/seasonality.py
@@ -28,11 +28,12 @@ class TimeSeasonality(Component):
         components are included in the same model. Default is ``f"Seasonal[s={season_length}, d={duration}]"``
 
     state_names: list of str, default None
-        List of strings for seasonal effect labels. If provided, it must be of length ``season_length``. An example
-        would be ``state_names = ['Mon', 'Tue', 'Wed', 'Thur', 'Fri', 'Sat', 'Sun']`` when data is daily with a weekly
+        List of strings for seasonal effect labels. If provided, it must be of length ``season_length`` times ``duration``.
+        An example would be ``state_names = ['Mon', 'Tue', 'Wed', 'Thur', 'Fri', 'Sat', 'Sun']`` when data is daily with a weekly
         seasonal pattern (``season_length = 7``).
 
-        If None, states will be numbered ``[State_0, ..., State_s-1]``
+        If None and ``duration = 1``, states will be named as ``[State_0, ..., State_s-1]`` (here s is ``season_length``).
+        If None and ``duration > 1``, states will be named as ``[State_0_0, ..., State_s-1_d-1]`` (here d is ``duration``).
 
     remove_first_state: bool, default True
         If True, the first state will be removed from the model. This is done because there are only ``season_length-1`` degrees of
@@ -148,7 +149,7 @@ class TimeSeasonality(Component):
     And so on. So for interpretation, the ``season_length - 1`` initial states are, when reversed, the coefficients
     associated with ``state_names[1:]``.
 
-    In the next example, we set :math:`s=3`, :math:`d=2`, ``remove_first_state=True``, and omit the shock term.
+    In the next example, we set :math:`s=2`, :math:`d=2`, ``remove_first_state=True``, and omit the shock term.
     By definition, the initial vector :math:`\alpha_{0}` is
 
     .. math::
@@ -178,7 +179,7 @@ class TimeSeasonality(Component):
 
     .. warning::
         Although the ``state_names`` argument expects a list of length ``season_length`` times ``duration``,
-        only ``state_names[1:]`` will be saved as model dimensions, since the first coefficient is not identified
+        only ``state_names[duration:]`` will be saved as model dimensions, since the first coefficient is not identified
         (it is defined as :math:`-\sum_{i=1}^{s-1} \tilde{\gamma}_{-i}`).
 
     Examples
@@ -269,7 +270,7 @@ def __init__(
 
         self.provided_state_names = state_names
 
-        k_states = season_length * duration - int(self.remove_first_state) * duration
+        k_states = (season_length - int(self.remove_first_state)) * duration
         k_endog = len(observed_state_names)
         k_posdef = int(innovations)
 

From 61373fe88f7c68d30b3163cf666aa58a1a0f86d8 Mon Sep 17 00:00:00 2001
From: gdeleva <giacomo.deleva@gmail.com>
Date: Fri, 18 Jul 2025 12:13:16 +0900
Subject: [PATCH 11/16] Handle initial states

---
 .../models/structural/components/seasonality.py     | 13 +++++++++++--
 1 file changed, 11 insertions(+), 2 deletions(-)

diff --git a/pymc_extras/statespace/models/structural/components/seasonality.py b/pymc_extras/statespace/models/structural/components/seasonality.py
index 9769ff95a..393a45ad8 100644
--- a/pymc_extras/statespace/models/structural/components/seasonality.py
+++ b/pymc_extras/statespace/models/structural/components/seasonality.py
@@ -332,6 +332,8 @@ def populate_component_properties(self):
 
     def make_symbolic_graph(self) -> None:
         k_states = self.k_states // self.k_endog
+        duration = self.duration
+        k_unique_states = k_states // duration
         k_posdef = self.k_posdef // self.k_endog
         k_endog = self.k_endog
 
@@ -364,6 +366,7 @@ def make_symbolic_graph(self) -> None:
             # circulant matrix that cycles between the states.
             T = np.eye(k_states, k=1)
             T[-1, 0] = 1
+        T = pt.as_tensor_variable(T)
 
         self.ssm["transition", :, :] = pt.linalg.block_diag(*[T for _ in range(k_endog)])
 
@@ -371,9 +374,15 @@ def make_symbolic_graph(self) -> None:
         self.ssm["design", :, :] = pt.linalg.block_diag(*[Z for _ in range(k_endog)])
 
         initial_states = self.make_and_register_variable(
-            f"coefs_{self.name}", shape=(k_states,) if k_endog == 1 else (k_endog, k_states)
+            f"coefs_{self.name}",
+            shape=(k_unique_states,) if k_endog == 1 else (k_endog, k_unique_states),
         )
-        self.ssm["initial_state", :] = initial_states.ravel()
+        if k_endog == 1:
+            self.ssm["initial_state", :] = pt.extra_ops.repeat(initial_states, duration, axis=0)
+        else:
+            self.ssm["initial_state", :] = pt.extra_ops.repeat(
+                initial_states, duration, axis=1
+            ).ravel()
 
         if self.innovations:
             R = pt.zeros((k_states, k_posdef))[0, 0].set(1.0)

From dc44ac5bd84685b03a873f891bea45e58dca7025 Mon Sep 17 00:00:00 2001
From: gdeleva <giacomo.deleva@gmail.com>
Date: Fri, 18 Jul 2025 22:48:56 +0900
Subject: [PATCH 12/16] Build symbolic Graph from pytensor operations

---
 .../models/structural/components/seasonality.py | 17 ++++++++---------
 1 file changed, 8 insertions(+), 9 deletions(-)

diff --git a/pymc_extras/statespace/models/structural/components/seasonality.py b/pymc_extras/statespace/models/structural/components/seasonality.py
index 393a45ad8..83027ec72 100644
--- a/pymc_extras/statespace/models/structural/components/seasonality.py
+++ b/pymc_extras/statespace/models/structural/components/seasonality.py
@@ -340,14 +340,14 @@ def make_symbolic_graph(self) -> None:
         if self.remove_first_state:
             # In this case, parameters are normalized to sum to zero, so the current state is the negative sum of
             # all previous states.
-            zero_d = np.zeros((self.duration, self.duration))
-            id_d = np.eye(self.duration)
+            zero_d = pt.zeros((self.duration, self.duration))
+            id_d = pt.eye(self.duration)
 
-            blocks = []
+            row_blocks = []
 
             # First row: all -1_d blocks
             first_row = [-id_d for _ in range(self.season_length - 1)]
-            blocks.append(first_row)
+            row_blocks.append(pt.concatenate(first_row, axis=1))
 
             # Rows 2 to season_length-1: shifted identity blocks
             for i in range(self.season_length - 2):
@@ -357,16 +357,15 @@ def make_symbolic_graph(self) -> None:
                         row.append(id_d)
                     else:
                         row.append(zero_d)
-                blocks.append(row)
+                row_blocks.append(pt.concatenate(row, axis=1))
 
             # Stack blocks
-            T = np.block(blocks)
+            T = pt.concatenate(row_blocks, axis=0)
         else:
             # In this case we assume the user to be responsible for ensuring the states sum to zero, so T is just a
             # circulant matrix that cycles between the states.
-            T = np.eye(k_states, k=1)
-            T[-1, 0] = 1
-        T = pt.as_tensor_variable(T)
+            T = pt.eye(k_states, k=1)
+            T = pt.set_subtensor(T[-1, 0], 1)
 
         self.ssm["transition", :, :] = pt.linalg.block_diag(*[T for _ in range(k_endog)])
 

From 9456bb2c6b6983bd2b9be56485648030e5b0e2d3 Mon Sep 17 00:00:00 2001
From: gdeleva <giacomo.deleva@gmail.com>
Date: Fri, 18 Jul 2025 22:50:03 +0900
Subject: [PATCH 13/16] Extend tests for the case duration > 1

---
 .../structural/components/test_seasonality.py | 138 ++++++++++--------
 1 file changed, 77 insertions(+), 61 deletions(-)

diff --git a/tests/statespace/models/structural/components/test_seasonality.py b/tests/statespace/models/structural/components/test_seasonality.py
index 93dcb0a5a..972a9c8c3 100644
--- a/tests/statespace/models/structural/components/test_seasonality.py
+++ b/tests/statespace/models/structural/components/test_seasonality.py
@@ -14,6 +14,7 @@
 
 
 @pytest.mark.parametrize("s", [10, 25, 50])
+@pytest.mark.parametrize("d", [1, 2, 3])
 @pytest.mark.parametrize("innovations", [True, False])
 @pytest.mark.parametrize("remove_first_state", [True, False])
 @pytest.mark.filterwarnings(
@@ -21,55 +22,61 @@
     "ignore:overflow encountered in matmul:RuntimeWarning",
     "ignore:invalid value encountered in matmul:RuntimeWarning",
 )
-def test_time_seasonality(s, innovations, remove_first_state, rng):
+def test_time_seasonality(s, d, innovations, remove_first_state, rng):
     def random_word(rng):
         return "".join(rng.choice(list("abcdefghijklmnopqrstuvwxyz")) for _ in range(5))
 
-    state_names = [random_word(rng) for _ in range(s)]
+    state_names = [random_word(rng) for _ in range(s * d)]
     mod = st.TimeSeasonality(
         season_length=s,
+        duration=d,
         innovations=innovations,
         name="season",
         state_names=state_names,
         remove_first_state=remove_first_state,
     )
-    x0 = np.zeros(mod.k_states, dtype=config.floatX)
+    x0 = np.zeros(mod.k_states // mod.duration, dtype=config.floatX)
     x0[0] = 1
 
     params = {"coefs_season": x0}
     if innovations:
         params["sigma_season"] = 0.0
 
-    x, y = simulate_from_numpy_model(mod, rng, params)
+    x, y = simulate_from_numpy_model(mod, rng, params, steps=100 * mod.duration)
     y = y.ravel()
     if not innovations:
-        assert_pattern_repeats(y, s, atol=ATOL, rtol=RTOL)
+        assert_pattern_repeats(y, s * d, atol=ATOL, rtol=RTOL)
 
     # Check coords
     mod = mod.build(verbose=False)
     _assert_basic_coords_correct(mod)
-    test_slice = slice(1, None) if remove_first_state else slice(None)
+    test_slice = slice(d, None) if remove_first_state else slice(None)
     assert mod.coords["state_season"] == state_names[test_slice]
 
 
+@pytest.mark.parametrize("d", [1, 2, 3])
 @pytest.mark.parametrize(
     "remove_first_state", [True, False], ids=["remove_first_state", "keep_first_state"]
 )
-def test_time_seasonality_multiple_observed(rng, remove_first_state):
+def test_time_seasonality_multiple_observed(rng, d, remove_first_state):
     s = 3
-    state_names = [f"state_{i}" for i in range(s)]
+    state_names = [f"state_{i}_{j}" for i in range(s) for j in range(d)]
     mod = st.TimeSeasonality(
         season_length=s,
+        duration=d,
         innovations=True,
         name="season",
         state_names=state_names,
         observed_state_names=["data_1", "data_2"],
         remove_first_state=remove_first_state,
     )
-    x0 = np.zeros((mod.k_endog, mod.k_states // mod.k_endog), dtype=config.floatX)
+    x0 = np.zeros((mod.k_endog, mod.k_states // mod.k_endog // mod.duration), dtype=config.floatX)
 
     expected_states = [
-        f"state_{i}[data_{j}]" for j in range(1, 3) for i in range(int(remove_first_state), s)
+        f"state_{i}_{j}[data_{k}]"
+        for k in range(1, 3)
+        for i in range(int(remove_first_state), s)
+        for j in range(d)
     ]
     assert mod.state_names == expected_states
     assert mod.shock_names == ["season[data_1]", "season[data_2]"]
@@ -79,9 +86,9 @@ def test_time_seasonality_multiple_observed(rng, remove_first_state):
 
     params = {"coefs_season": x0, "sigma_season": np.array([0.0, 0.0], dtype=config.floatX)}
 
-    x, y = simulate_from_numpy_model(mod, rng, params, steps=123)
-    assert_pattern_repeats(y[:, 0], s, atol=ATOL, rtol=RTOL)
-    assert_pattern_repeats(y[:, 1], s, atol=ATOL, rtol=RTOL)
+    x, y = simulate_from_numpy_model(mod, rng, params, steps=123 * d)
+    assert_pattern_repeats(y[:, 0], s * d, atol=ATOL, rtol=RTOL)
+    assert_pattern_repeats(y[:, 1], s * d, atol=ATOL, rtol=RTOL)
 
     mod = mod.build(verbose=False)
     x0, *_, T, Z, R, _, Q = mod._unpack_statespace_with_placeholders()
@@ -97,36 +104,38 @@ def test_time_seasonality_multiple_observed(rng, remove_first_state):
     params["sigma_season"] = np.array([0.1, 0.8], dtype=config.floatX)
     x0, T, Z, R, Q = fn(**params)
 
+    # Because the dimension of the observed states is 2,
+    # the expected T is the diagonal block matrix [[T0, 0], [0, T0]]
+    # where T0 is the transition matrix we would have if the
+    # seasonality were not multiple observed.
+    mod0 = st.TimeSeasonality(season_length=s, duration=d, remove_first_state=remove_first_state)
+    T0 = mod0.ssm["transition"].eval()
+
     if remove_first_state:
-        expected_x0 = np.array([1.0, 0.0, 2.0, 0.0])
-
-        expected_T = np.array(
-            [
-                [-1.0, -1.0, 0.0, 0.0],
-                [1.0, 0.0, 0.0, 0.0],
-                [0.0, 0.0, -1.0, -1.0],
-                [0.0, 0.0, 1.0, 0.0],
-            ]
+        expected_x0 = np.repeat(np.array([1.0, 0.0, 2.0, 0.0]), d)
+        expected_T = np.block(
+            [[T0, np.zeros((d * (s - 1), d * (s - 1)))], [np.zeros((d * (s - 1), d * (s - 1))), T0]]
+        )
+        expected_R = np.array(
+            [[1.0, 1.0]] + [[0.0, 0.0]] * (2 * d - 1) + [[1.0, 1.0]] + [[0.0, 0.0]] * (2 * d - 1)
         )
-        expected_R = np.array([[1.0, 1.0], [0.0, 0.0], [1.0, 1.0], [0.0, 0.0]])
-        expected_Z = np.array([[1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0]])
+        Z0 = np.zeros((2, d * (s - 1)))
+        Z0[0, 0] = 1
+        Z1 = np.zeros((2, d * (s - 1)))
+        Z1[1, 0] = 1
+        expected_Z = np.block([[Z0, Z1]])
 
     else:
-        expected_x0 = np.array([1.0, 0.0, 0.0, 2.0, 0.0, 0.0])
-        expected_T = np.array(
-            [
-                [0.0, 1.0, 0.0, 0.0, 0.0, 0.0],
-                [0.0, 0.0, 1.0, 0.0, 0.0, 0.0],
-                [1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
-                [0.0, 0.0, 0.0, 0.0, 1.0, 0.0],
-                [0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
-                [0.0, 0.0, 0.0, 1.0, 0.0, 0.0],
-            ]
-        )
+        expected_x0 = np.repeat(np.array([1.0, 0.0, 0.0, 2.0, 0.0, 0.0]), d)
+        expected_T = np.block([[T0, np.zeros((s * d, s * d))], [np.zeros((s * d, s * d)), T0]])
         expected_R = np.array(
-            [[1.0, 1.0], [0.0, 0.0], [0.0, 0.0], [1.0, 1.0], [0.0, 0.0], [0.0, 0.0]]
+            [[1.0, 1.0]] + [[0.0, 0.0]] * (s * d - 1) + [[1.0, 1.0]] + [[0.0, 0.0]] * (s * d - 1)
         )
-        expected_Z = np.array([[1.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0, 0.0, 0.0]])
+        Z0 = np.zeros((2, s * d))
+        Z0[0, 0] = 1
+        Z1 = np.zeros((2, s * d))
+        Z1[1, 0] = 1
+        expected_Z = np.block([[Z0, Z1]])
 
     expected_Q = np.array([[0.1**2, 0.0], [0.0, 0.8**2]])
 
@@ -137,20 +146,24 @@ def test_time_seasonality_multiple_observed(rng, remove_first_state):
         np.testing.assert_allclose(matrix, expected)
 
 
-def test_add_two_time_seasonality_different_observed(rng):
+@pytest.mark.parametrize("d1", [1, 2, 3])
+@pytest.mark.parametrize("d2", [1, 2, 3])
+def test_add_two_time_seasonality_different_observed(rng, d1, d2):
     mod1 = st.TimeSeasonality(
         season_length=3,
+        duration=d1,
         innovations=True,
         name="season1",
-        state_names=[f"state_{i}" for i in range(3)],
+        state_names=[f"state_{i}_{j}" for i in range(3) for j in range(d1)],
         observed_state_names=["data_1"],
         remove_first_state=False,
     )
     mod2 = st.TimeSeasonality(
         season_length=5,
+        duration=d2,
         innovations=True,
         name="season2",
-        state_names=[f"state_{i}" for i in range(5)],
+        state_names=[f"state_{i}_{j}" for i in range(5) for j in range(d2)],
         observed_state_names=["data_2"],
     )
 
@@ -164,18 +177,22 @@ def test_add_two_time_seasonality_different_observed(rng):
         "initial_state_cov": np.eye(mod.k_states, dtype=config.floatX),
     }
 
-    x, y = simulate_from_numpy_model(mod, rng, params, steps=3 * 5 * 5)
-    assert_pattern_repeats(y[:, 0], 3, atol=ATOL, rtol=RTOL)
-    assert_pattern_repeats(y[:, 1], 5, atol=ATOL, rtol=RTOL)
+    x, y = simulate_from_numpy_model(mod, rng, params, steps=3 * 5 * 5 * d1 * d2)
+    assert_pattern_repeats(y[:, 0], 3 * d1, atol=ATOL, rtol=RTOL)
+    assert_pattern_repeats(y[:, 1], 5 * d2, atol=ATOL, rtol=RTOL)
 
     assert mod.state_names == [
-        "state_0[data_1]",
-        "state_1[data_1]",
-        "state_2[data_1]",
-        "state_1[data_2]",
-        "state_2[data_2]",
-        "state_3[data_2]",
-        "state_4[data_2]",
+        item
+        for sublist in [
+            [f"state_0_{j}[data_1]" for j in range(d1)],
+            [f"state_1_{j}[data_1]" for j in range(d1)],
+            [f"state_2_{j}[data_1]" for j in range(d1)],
+            [f"state_1_{j}[data_2]" for j in range(d2)],
+            [f"state_2_{j}[data_2]" for j in range(d2)],
+            [f"state_3_{j}[data_2]" for j in range(d2)],
+            [f"state_4_{j}[data_2]" for j in range(d2)],
+        ]
+        for item in sublist
     ]
 
     assert mod.shock_names == ["season1[data_1]", "season2[data_2]"]
@@ -194,20 +211,19 @@ def test_add_two_time_seasonality_different_observed(rng):
     )
 
     np.testing.assert_allclose(
-        np.array([1.0, 0.0, 0.0, 3.0, 0.0, 0.0, 1.2]), x0, atol=ATOL, rtol=RTOL
+        np.repeat(np.array([1.0, 0.0, 0.0, 3.0, 0.0, 0.0, 1.2]), [d1, d1, d1, d2, d2, d2, d2]),
+        x0,
+        atol=ATOL,
+        rtol=RTOL,
     )
 
+    # The transition matrix T of mod is expected to be [[T1, 0], [0, T2]],
+    # where T1 and T2 are the transition matrices of mod1 and mod2, respectively.
+    T1 = mod1.ssm["transition"].eval()
+    T2 = mod2.ssm["transition"].eval()
     np.testing.assert_allclose(
-        np.array(
-            [
-                [0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
-                [0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0],
-                [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
-                [0.0, 0.0, 0.0, -1.0, -1.0, -1.0, -1.0],
-                [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0],
-                [0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0],
-                [0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0],
-            ]
+        np.block(
+            [[T1, np.zeros((T1.shape[0], T2.shape[1]))], [np.zeros((T2.shape[0], T1.shape[1])), T2]]
         ),
         T,
         atol=ATOL,

From 02490e78ca3fe4611ac212d87d40cd0cb4983971 Mon Sep 17 00:00:00 2001
From: gdeleva <giacomo.deleva@gmail.com>
Date: Mon, 21 Jul 2025 20:40:05 +0900
Subject: [PATCH 14/16] Reduce number of tests for duration

---
 .../models/structural/components/test_seasonality.py       | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/tests/statespace/models/structural/components/test_seasonality.py b/tests/statespace/models/structural/components/test_seasonality.py
index 972a9c8c3..162b25868 100644
--- a/tests/statespace/models/structural/components/test_seasonality.py
+++ b/tests/statespace/models/structural/components/test_seasonality.py
@@ -14,7 +14,7 @@
 
 
 @pytest.mark.parametrize("s", [10, 25, 50])
-@pytest.mark.parametrize("d", [1, 2, 3])
+@pytest.mark.parametrize("d", [1, 3])
 @pytest.mark.parametrize("innovations", [True, False])
 @pytest.mark.parametrize("remove_first_state", [True, False])
 @pytest.mark.filterwarnings(
@@ -54,7 +54,7 @@ def random_word(rng):
     assert mod.coords["state_season"] == state_names[test_slice]
 
 
-@pytest.mark.parametrize("d", [1, 2, 3])
+@pytest.mark.parametrize("d", [1, 3])
 @pytest.mark.parametrize(
     "remove_first_state", [True, False], ids=["remove_first_state", "keep_first_state"]
 )
@@ -146,8 +146,7 @@ def test_time_seasonality_multiple_observed(rng, d, remove_first_state):
         np.testing.assert_allclose(matrix, expected)
 
 
-@pytest.mark.parametrize("d1", [1, 2, 3])
-@pytest.mark.parametrize("d2", [1, 2, 3])
+@pytest.mark.parametrize("d1, d2", [(1, 1), (1, 3), (3, 1), (3, 3)])
 def test_add_two_time_seasonality_different_observed(rng, d1, d2):
     mod1 = st.TimeSeasonality(
         season_length=3,

From c3e74fbe513cb82d5f403b059bafd206ae5d46fa Mon Sep 17 00:00:00 2001
From: gdeleva <giacomo.deleva@gmail.com>
Date: Mon, 21 Jul 2025 21:49:50 +0900
Subject: [PATCH 15/16] Merge & add conditions for raising ValueErrors for
 season_length and duration

---
 .../models/structural/components/seasonality.py           | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/pymc_extras/statespace/models/structural/components/seasonality.py b/pymc_extras/statespace/models/structural/components/seasonality.py
index 83027ec72..4b9d776d4 100644
--- a/pymc_extras/statespace/models/structural/components/seasonality.py
+++ b/pymc_extras/statespace/models/structural/components/seasonality.py
@@ -239,8 +239,12 @@ def __init__(
         if observed_state_names is None:
             observed_state_names = ["data"]
 
-        if season_length <= 1:
-            raise ValueError(f"season_length must be greater than 1, got {season_length}")
+        if season_length <= 1 or not isinstance(duration, int):
+            raise ValueError(
+                f"season_length must be an integer greater than 1, got {season_length}"
+            )
+        if duration <= 0 or not isinstance(duration, int):
+            raise ValueError(f"duration must be a positive integer, got {duration}")
         if name is None:
             name = f"Seasonal[s={season_length}, d={duration}]"
         if state_names is None:

From 3bb28fd166038e32b81ba5874ef10ae2d955259f Mon Sep 17 00:00:00 2001
From: gdeleva <giacomo.deleva@gmail.com>
Date: Mon, 21 Jul 2025 21:51:17 +0900
Subject: [PATCH 16/16] Fix bug

---
 .../statespace/models/structural/components/seasonality.py      | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pymc_extras/statespace/models/structural/components/seasonality.py b/pymc_extras/statespace/models/structural/components/seasonality.py
index 4b9d776d4..341bd6ae5 100644
--- a/pymc_extras/statespace/models/structural/components/seasonality.py
+++ b/pymc_extras/statespace/models/structural/components/seasonality.py
@@ -239,7 +239,7 @@ def __init__(
         if observed_state_names is None:
             observed_state_names = ["data"]
 
-        if season_length <= 1 or not isinstance(duration, int):
+        if season_length <= 1 or not isinstance(season_length, int):
             raise ValueError(
                 f"season_length must be an integer greater than 1, got {season_length}"
             )