Feature/improved cyclical (#30)

* Adds reference

* Adds cyclical

* Fixes cyclical, error in statsmodels docs

* Adds test

* Meta

CHANGELOG.md

@@ -2,6 +2,9 @@
 
 ## Versions
 
+### v0.2.4
+ - Adds `CyclicalProcess` and "fixes" `HarmonicProcess`
+
 ### v0.2.3
  - Bug fixes for `DoubleExponential` and `SelfExcitingProcess`.
 

pyproject.toml

@@ -46,7 +46,7 @@ line-length = 120
 include = ["stochproc*"]
 
 [tool.bumpver]
-current_version = "0.2.3"
+current_version = "0.2.4"
 version_pattern = "MAJOR.MINOR.PATCH"
 commit_message = "bump version {old_version} -> {new_version}"
 commit = false

stochproc/__init__.py

@@ -1,4 +1,4 @@
-__version__ = "0.2.3"
+__version__ = "0.2.4"
 
 from . import distributions
 from . import timeseries

stochproc/timeseries/models/__init__.py

@@ -8,8 +8,10 @@
 from .smooth_trend import SmoothLinearTrend
 from .trending_ou import TrendingOU
 from .self_exciting_process import SelfExcitingLatentProcesses
+from .cyclical import CyclicalProcess
 from .harmonics import HarmonicProcess
 
+
 __all__ = [
     "OrnsteinUhlenbeck",
     "Verhulst",
@@ -22,4 +24,5 @@
     "TrendingOU",
     "SelfExcitingLatentProcesses",
     "HarmonicProcess",
+    "CyclicalProcess",
 ]

stochproc/timeseries/models/cyclical.py

@@ -0,0 +1,52 @@
+import torch
+from pyro.distributions import Delta, Normal
+from torch.distributions.utils import broadcast_all
+
+from ...typing import ParameterType
+from ..linear import LinearModel
+
+
+def _initial_kernel(v_0):
+    return Delta(v_0, event_dim=1)
+
+
+def _parameter_transform(rho, lamda, s):
+    cos_lam = rho * torch.cos(lamda)
+    sin_lam = rho * torch.sin(lamda)
+
+    a_top = torch.stack([cos_lam, sin_lam], dim=-1)
+    a_bottom = torch.stack([-sin_lam, cos_lam], dim=-1)
+
+    a = torch.stack([a_top, a_bottom], dim=-2)
+
+    return a, torch.zeros_like(s), s
+
+
+class CyclicalProcess(LinearModel):
+    """
+    Implements a cyclical process like `statsmodels`_.
+
+    .. _`statsmodels`: https://www.statsmodels.org/stable/generated/statsmodels.tsa.statespace.structural.UnobservedComponents.html#statsmodels.tsa.statespace.structural.UnobservedComponents
+    """
+
+    def __init__(self, rho: ParameterType, lamda: ParameterType, sigma: ParameterType, x_0: ParameterType = None):
+        """
+        Internal initializer for :class:`Cyclical`.
+
+        Args:
+            rho (ParameterType): see reference.
+            lamda (ParameterType): see reference.
+            sigma (ParameterType): see reference.
+            x_0 (ParameterType): initial values.
+        """
+
+        rho, lamda, sigma = broadcast_all(rho, lamda, sigma)
+
+        if x_0 is None:
+            x_0 = torch.zeros(2, device=lamda.device)
+
+        distribution = Normal(
+            torch.tensor(0.0, device=rho.device), torch.tensor(1.0, device=rho.device)
+            ).expand(torch.Size([2])).to_event(1)
+
+        super().__init__((rho, lamda, sigma), distribution, _initial_kernel, initial_parameters=(x_0,), parameter_transform=_parameter_transform)

stochproc/timeseries/models/harmonics.py

@@ -1,60 +1,32 @@
-from pyro.distributions import Normal
-import torch
+from math import pi
 
-from ..linear import LinearModel
+from .cyclical import CyclicalProcess
 from ..utils import coerce_tensors
 from ...typing import ParameterType
 
 
-def _parameter_transform(lamda, s):
-    cos_lam = torch.cos(lamda)
-    sin_lam = torch.sin(lamda)
-
-    a_top = torch.stack([cos_lam, sin_lam], dim=-1)
-    a_bottom = torch.stack([-sin_lam, cos_lam], dim=-1)
-
-    a = torch.stack([a_top, a_bottom], dim=-2)
-
-    return a, torch.zeros_like(s), s
-
-
-def initial_kernel(x0, s):
-    return Normal(x0, s).to_event(1)
-
-
-class HarmonicProcess(LinearModel):
+class HarmonicProcess(CyclicalProcess):
     r"""
     Implements a harmonic timeseries process of the form
         .. math::
             \gamma_{t + 1} = \gamma \cos{ \lambda } + \gamma^*\sin{ \lambda } + \sigma \nu_{t + 1}, \newline
             \gamma^*_{t + 1} = -\gamma \sin { \lambda } + \gamma^* \cos{ \lambda } + \sigma^* \nu^*_{t + 1}.
+        
+    See `statsmodels`_.
+
+    .. _`statsmodels`: https://www.statsmodels.org/stable/generated/statsmodels.tsa.statespace.structural.UnobservedComponents.html#statsmodels.tsa.statespace.structural.UnobservedComponents
     """
 
-    def __init__(self, lamda: ParameterType, sigma: ParameterType, x_0: ParameterType = None):
+    def __init__(self, s: int, sigma: ParameterType, x_0: ParameterType = None, j: int = 1):
         """
         Internal initializer for :class:`HarmonicProcess`.
 
         Args:
-            lamda (ParameterType): coefficient for periodic component.
-            sigma (ParameterType): the st
+            s (int): number of seasons.
+            sigma (ParameterType): see :class:`Cyclical`.
+            x_0 (ParameterType): see :class:`Cyclical`.
+            j (int): "index" of harmonic process.
         """
 
-        lamda, sigma = coerce_tensors(lamda, sigma)
-
-        if x_0 is None:
-            x_0 = torch.zeros(2, device=lamda.device)
-
-        lamda, sigma, x_0 = coerce_tensors(lamda, sigma, x_0)
-        increment_distribution = Normal(
-            torch.zeros(2, device=lamda.device), torch.ones(2, device=lamda.device)
-        ).to_event(1)
-
-        initial_parameters = (x_0, sigma)
-
-        super().__init__(
-            (lamda, sigma),
-            increment_distribution,
-            initial_kernel,
-            initial_parameters=initial_parameters,
-            parameter_transform=_parameter_transform,
-        )
+        rho, lamda, sigma = coerce_tensors(1.0, 2.0 * pi * j / s, sigma)
+        super().__init__(rho, lamda, sigma, x_0)

tests/timeseries/test_custom.py

@@ -1,4 +1,5 @@
 import itertools
+from math import pi
 
 import pytest
 import torch
@@ -21,7 +22,8 @@ def models():
     yield mods.SmoothLinearTrend(mods.OrnsteinUhlenbeck(0.025, 0.0, 0.05))
     yield mods.TrendingOU(0.01, 0.03, 0.2, 1.0)
     yield mods.SelfExcitingLatentProcesses(0.01, 2.0, 0.05, 0.1, 3.0, 2.0, dt=0.05)
-    yield mods.HarmonicProcess(0.25, torch.tensor([0.05, 0.05]))
+    yield mods.HarmonicProcess(3, 0.05)
+    yield mods.CyclicalProcess(0.98, 2.0 * pi / 1_000, 0.25)
 
 
 class TestCustomModels(object):