Hamilton filter updates

statsmodels/tsa/regime_switching/_hamilton_filter.pyx.in

@@ -4,7 +4,7 @@
 """
 Hamilton filter
 
-Author: Chad Fulton  
+Author: Chad Fulton
 License: Simplified-BSD
 """
 
@@ -38,29 +38,32 @@ def {{prefix}}hamilton_filter(int nobs, int k_regimes, int order,
                               {{cython_type}} [:,:] filtered_joint_probabilities):
     cdef int t, i, j, k, ix, regime_transition_t = 0, time_varying_regime_transition
     cdef:
-        int k_regimes_order_m1 = k_regimes**(order - 1)
+        # k_regimes_order_m1 is not used when order == 0.
+        int k_regimes_order_m1 = k_regimes**max(order - 1, 0)
         int k_regimes_order = k_regimes**order
         int k_regimes_order_p1 = k_regimes**(order + 1)
         {{cython_type}} [:] weighted_likelihoods, tmp_filtered_marginalized_probabilities
 
     time_varying_regime_transition = regime_transition.shape[2] > 1
     weighted_likelihoods = np.zeros(k_regimes_order_p1, dtype={{dtype}})
+    # tmp_filtered_marginalized_probabilities is not used if order == 0.
     tmp_filtered_marginalized_probabilities = np.zeros(k_regimes_order, dtype={{dtype}})
 
     for t in range(nobs):
         if time_varying_regime_transition:
             regime_transition_t = t
 
-        # Collapse filtered joint probabilities over the last dimension
-        # Pr[S_{t-1}, ..., S_{t-r} | t-1] = \sum_{ S_{t-r-1} } Pr[S_{t-1}, ..., S_{t-r}, S_{t-r-1} | t-1]
-        ix = 0
-        tmp_filtered_marginalized_probabilities[:] = 0
-        for j in range(k_regimes_order):
-            for i in range(k_regimes):
-                tmp_filtered_marginalized_probabilities[j] = (
-                    tmp_filtered_marginalized_probabilities[j] +
-                    filtered_joint_probabilities[ix, t])
-                ix = ix + 1
+        if order > 0:
+            # Collapse filtered joint probabilities over the last dimension
+            # Pr[S_{t-1}, ..., S_{t-r} | t-1] = \sum_{ S_{t-r-1} } Pr[S_{t-1}, ..., S_{t-r}, S_{t-r-1} | t-1]
+            ix = 0
+            tmp_filtered_marginalized_probabilities[:] = 0
+            for j in range(k_regimes_order):
+                for i in range(k_regimes):
+                    tmp_filtered_marginalized_probabilities[j] = (
+                        tmp_filtered_marginalized_probabilities[j] +
+                        filtered_joint_probabilities[ix, t])
+                    ix = ix + 1
 
         {{prefix}}hamilton_filter_iteration(t, k_regimes, order,
                                   regime_transition[:, :, regime_transition_t],
@@ -90,14 +93,25 @@ cdef {{prefix}}hamilton_filter_iteration(int t, int k_regimes, int order,
 
     # Compute predicted joint probabilities
     # Pr[S_t, S_{t-1}, ..., S_{t-r} | t-1] = Pr[S_t | S_{t-1}] * Pr[S_{t-1}, ..., S_{t-r} | t-1]
-    ix = 0
-    for i in range(k_regimes):
-        for j in range(k_regimes):
-            for k in range(k_regimes_order_m1):
-                curr_predicted_joint_probabilities[ix] = (
-                    prev_filtered_marginalized_probabilities[j * k_regimes_order_m1 + k] *
-                    regime_transition[i, j])
-                ix += 1
+    if order > 0:
+        ix = 0
+        for i in range(k_regimes):
+            for j in range(k_regimes):
+                for k in range(k_regimes_order_m1):
+                    curr_predicted_joint_probabilities[ix] = (
+                        prev_filtered_marginalized_probabilities[j * k_regimes_order_m1 + k] *
+                        regime_transition[i, j])
+                    ix += 1
+    else:
+        curr_predicted_joint_probabilities[:] = 0
+        for i in range(k_regimes):
+            for j in range(k_regimes):
+                # There appears to be a bug in cython for += with complex types.
+                # https://groups.google.com/forum/#!topic/cython-users/jD8U6AuYKS0
+                curr_predicted_joint_probabilities[i] = (
+                    curr_predicted_joint_probabilities[i]
+                    + prev_filtered_joint_probabilities[j] * regime_transition[i, j])
+
 
     # Compute weighted likelihoods f(y_t | S_t, S_{t-1}, ..., S_{t-r}, t-1) * Pr[S_t, S_{t-1}, ..., S_{t-r} | t-1]
     # and the joint likelihood f(y_t | t-1)

statsmodels/tsa/regime_switching/markov_switching.py

@@ -110,11 +110,18 @@ def py_hamilton_filter(initial_probabilities, regime_transition,
     Parameters
     ----------
     initial_probabilities : array
-        Array of initial probabilities, shaped (k_regimes,).
+        Array of initial probabilities, shaped (k_regimes,) giving the
+        distribution of the regime process at time t = -order where order
+        is a nonnegative integer.
     regime_transition : array
         Matrix of regime transition probabilities, shaped either
         (k_regimes, k_regimes, 1) or if there are time-varying transition
-        probabilities (k_regimes, k_regimes, nobs).
+        probabilities (k_regimes, k_regimes, nobs + order).  Entry [i, j,
+        t] contains the probability of moving from j at time t-1 to i at
+        time t, so each matrix regime_transition[:, :, t] should be left
+        stochastic.  The first order entries and initial_probabilities are
+        used to produce the initial joint distribution of dimension order +
+        1 at time t=0.
     conditional_likelihoods : array
         Array of likelihoods conditional on the last `order+1` regimes,
         shaped (k_regimes,)*(order + 1) + (nobs,).
@@ -136,6 +143,7 @@ def py_hamilton_filter(initial_probabilities, regime_transition,
         the joint probability of the current and previous `order` periods
         being in each combination of regimes conditional on time t
         information. Shaped (k_regimes,) * (order + 1) + (nobs,).
+
     """
 
     # Dimensions
@@ -144,6 +152,14 @@ def py_hamilton_filter(initial_probabilities, regime_transition,
     order = conditional_likelihoods.ndim - 2
     dtype = conditional_likelihoods.dtype
 
+    # Check for compatible shapes.
+    incompatible_shapes = (
+        regime_transition.shape[-1] not in (1, nobs + order)
+        or regime_transition.shape[:2] != (k_regimes, k_regimes)
+        or conditional_likelihoods.shape[0] != k_regimes)
+    if incompatible_shapes:
+        raise ValueError('Arguments do not have compatible shapes')
+
     # Storage
     # Pr[S_t = s_t | Y_t]
     filtered_marginal_probabilities = (
@@ -165,6 +181,11 @@ def py_hamilton_filter(initial_probabilities, regime_transition,
         tmp = np.reshape(regime_transition[..., i], shape + (1,) * i) * tmp
     filtered_joint_probabilities[..., 0] = tmp
 
+    # Check that regime_transition is oriented correctly.
+    if not np.allclose(np.sum(regime_transition, axis=0), 1):
+        raise ValueError('regime_transition does not contain '
+                         'left stochastic matrices.')
+
     # Reshape regime_transition so we can use broadcasting
     shape = (k_regimes, k_regimes)
     shape += (1,) * (order-1)
@@ -182,11 +203,16 @@ def py_hamilton_filter(initial_probabilities, regime_transition,
             transition_t = t
 
         # S_t, S_{t-1}, ..., S_{t-r} | t-1, stored at zero-indexed location t
-        predicted_joint_probabilities[..., t] = (
-            # S_t | S_{t-1}
-            regime_transition[..., transition_t] *
-            # S_{t-1}, S_{t-2}, ..., S_{t-r} | t-1
-            filtered_joint_probabilities[..., t].sum(axis=-1))
+        if order > 0:
+            predicted_joint_probabilities[..., t] = (
+                # S_t | S_{t-1}
+                regime_transition[..., transition_t] *
+                # S_{t-1}, S_{t-2}, ..., S_{t-r} | t-1
+                filtered_joint_probabilities[..., t].sum(axis=-1))
+        else:
+            predicted_joint_probabilities[..., t] = (
+                np.dot(regime_transition[..., transition_t],
+                       filtered_joint_probabilities[..., t]))
 
         # f(y_t, S_t, ..., S_{t-r} | t-1)
         tmp = (conditional_likelihoods[..., t] *
@@ -216,11 +242,18 @@ def cy_hamilton_filter(initial_probabilities, regime_transition,
     Parameters
     ----------
     initial_probabilities : array
-        Array of initial probabilities, shaped (k_regimes,).
+        Array of initial probabilities, shaped (k_regimes,) giving the
+        distribution of the regime process at time t = -order where order
+        is a nonnegative integer.
     regime_transition : array
         Matrix of regime transition probabilities, shaped either
         (k_regimes, k_regimes, 1) or if there are time-varying transition
-        probabilities (k_regimes, k_regimes, nobs).
+        probabilities (k_regimes, k_regimes, nobs + order).  Entry [i, j,
+        t] contains the probability of moving from j at time t-1 to i at
+        time t, so each matrix regime_transition[:, :, t] should be left
+        stochastic.  The first order entries and initial_probabilities are
+        used to produce the initial joint distribution of dimension order +
+        1 at time t=0.
     conditional_likelihoods : array
         Array of likelihoods conditional on the last `order+1` regimes,
         shaped (k_regimes,)*(order + 1) + (nobs,).
@@ -250,6 +283,14 @@ def cy_hamilton_filter(initial_probabilities, regime_transition,
     order = conditional_likelihoods.ndim - 2
     dtype = conditional_likelihoods.dtype
 
+    # Check for compatible shapes.
+    incompatible_shapes = (
+        regime_transition.shape[-1] not in (1, nobs + order)
+        or regime_transition.shape[:2] != (k_regimes, k_regimes)
+        or conditional_likelihoods.shape[0] != k_regimes)
+    if incompatible_shapes:
+        raise ValueError('Arguments do not have compatible shapes')
+
     # Storage
     # Pr[S_t = s_t | Y_t]
     filtered_marginal_probabilities = (

statsmodels/tsa/regime_switching/tests/test_markov_regression.py

@@ -12,7 +12,7 @@
 import pandas as pd
 from statsmodels.tsa.regime_switching import (markov_switching,
                                               markov_regression)
-from numpy.testing import assert_allclose
+from numpy.testing import assert_allclose, assert_raises
 
 
 current_path = os.path.dirname(os.path.abspath(__file__))
@@ -939,6 +939,102 @@ def test_smoother_output(self, **kwargs):
         assert_allclose(res.smoothed_joint_probabilities,
                         fedfunds_const_short_smoothed_joint_probabilities)
 
+    def test_hamilton_filter_order_zero(self):
+        k_regimes = 3
+        nobs = 4
+        initial_probabilities = np.ones(k_regimes) / k_regimes
+
+        # We don't actually transition between the 3 regimes.
+        regime_transition = np.eye(k_regimes)[:, :, np.newaxis]
+
+        # Regime i correponds to a sequence of iid draws from discrete
+        # random variables that are equally likely to be -1 and i for i=0,
+        # 1, 2.  We observe the sequence -1, -1, 1, -1.  The first two
+        # observations tell us nothing, but the third lets us know that we
+        # are in regime 1 the whole time.
+        conditional_likelihoods = np.ones((k_regimes, nobs)) / 2
+        conditional_likelihoods[:, 2] = [0, 1, 0]
+
+        expected_marginals = np.empty((k_regimes, nobs))
+        expected_marginals[:, :2] = [[1/3], [1/3], [1/3]]
+        expected_marginals[:, 2:] = [[0], [1], [0]]
+
+
+        py_results = markov_switching.py_hamilton_filter(
+            initial_probabilities, regime_transition, conditional_likelihoods)
+        assert_allclose(py_results[0], expected_marginals)
+
+        cy_results = markov_switching.cy_hamilton_filter(
+            initial_probabilities, regime_transition, conditional_likelihoods)
+        assert_allclose(cy_results[0], expected_marginals)
+
+    def test_hamilton_filter_order_zero_with_tvtp(self):
+        k_regimes = 3
+        nobs = 8
+        initial_probabilities = np.ones(k_regimes) / k_regimes
+
+        # We don't actually transition between the 3 regimes except from
+        # t=3 to t=4 where we reset to regimes 1 and 2 being equally
+        # likely.
+        regime_transition = np.zeros((k_regimes, k_regimes, nobs))
+        regime_transition[...] = np.eye(k_regimes)[:, :, np.newaxis]
+        regime_transition[..., 4] = [[0,     0,   0],
+                                     [1/2, 1/2, 1/2],
+                                     [1/2, 1/2, 1/2]]
+
+        # Regime i correponds to a sequence of iid draws from discrete
+        # random variables that are equally likely to be -1, i, and i + 1
+        # for i = 0, 1, 2.  We observe the sequence:
+        #
+        #     t =  0, 1, 2,  3,  4, 5, 6,  7
+        #   X_t = -1, 1, 2, -1, -1, 2, 3, -1
+        #
+        # The first observations tell us nothing, the second tells us that
+        # regime 0 and 1 are equally likely, and the third tells us that we
+        # must be in regime 1 for t = 0, 1, 2, 3.  At t=4 we transition to
+        # state 1 or 2.  In this case, neither a -1 or 2 changes our view
+        # that these states are equally likely, but a 3 tells we must be in
+        # state 2 for the second four timestamps.
+        conditional_likelihoods = np.empty((k_regimes, nobs))
+        conditional_likelihoods[:, 0] = [1/3, 1/3, 1/3]
+        conditional_likelihoods[:, 1] = [1/3, 1/3, 0]
+        conditional_likelihoods[:, 2] = [0, 1/3, 1/3]
+        conditional_likelihoods[:, 3:5] = [[1/3], [1/3], [1/3]]
+        conditional_likelihoods[:, 5] = [0, 1/3, 1/3]
+        conditional_likelihoods[:, 6] = [0, 0, 1/3]
+        conditional_likelihoods[:, 7] = [1/3, 1/3, 1/3]
+
+        expected_marginals = np.empty((k_regimes, nobs))
+        expected_marginals[:, 0] = [1/3, 1/3, 1/3]
+        expected_marginals[:, 1] = [1/2, 1/2, 0]
+        expected_marginals[:, 2:4] = [[0], [1], [0]]
+        expected_marginals[:, 4:6] = [[0], [1/2], [1/2]]
+        expected_marginals[:, 6:8] = [[0], [0], [1]]
+
+        py_results = markov_switching.py_hamilton_filter(
+            initial_probabilities, regime_transition, conditional_likelihoods)
+        assert_allclose(py_results[0], expected_marginals)
+
+        cy_results = markov_switching.cy_hamilton_filter(
+            initial_probabilities, regime_transition, conditional_likelihoods)
+        assert_allclose(cy_results[0], expected_marginals)
+
+    def test_hamilton_filter_shape_checks(self):
+        k_regimes = 3
+        nobs = 8
+        order = 3
+
+        initial_probabilities = np.ones(k_regimes) / k_regimes
+        regime_transition = np.ones((k_regimes, k_regimes, nobs)) / k_regimes
+        conditional_likelihoods = np.ones(order * (k_regimes,) + (nobs,))
+
+        for func in [markov_switching.py_hamilton_filter,
+                     markov_switching.cy_hamilton_filter]:
+            with assert_raises(ValueError):
+                func(initial_probabilities,
+                     regime_transition,
+                     conditional_likelihoods)
+
     def test_py_hamilton_filter(self):
         mod = self.model
         params = self.true['params']