From 0f06b1078e5ec634e2e979d9414a61569b6f2ba5 Mon Sep 17 00:00:00 2001 From: Gerard Brunick Date: Tue, 27 Feb 2018 20:33:58 -0500 Subject: [PATCH 1/3] Handle order=0 case in Hamilton filter --- .../regime_switching/_hamilton_filter.pyx.in | 54 ++++++++----- .../tsa/regime_switching/markov_switching.py | 15 ++-- .../tests/test_markov_regression.py | 80 +++++++++++++++++++ 3 files changed, 124 insertions(+), 25 deletions(-) diff --git a/statsmodels/tsa/regime_switching/_hamilton_filter.pyx.in b/statsmodels/tsa/regime_switching/_hamilton_filter.pyx.in index 73a21235590..6bb302fd7da 100644 --- a/statsmodels/tsa/regime_switching/_hamilton_filter.pyx.in +++ b/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) diff --git a/statsmodels/tsa/regime_switching/markov_switching.py b/statsmodels/tsa/regime_switching/markov_switching.py index 1a55330ee90..9182faa4d40 100644 --- a/statsmodels/tsa/regime_switching/markov_switching.py +++ b/statsmodels/tsa/regime_switching/markov_switching.py @@ -182,11 +182,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] * diff --git a/statsmodels/tsa/regime_switching/tests/test_markov_regression.py b/statsmodels/tsa/regime_switching/tests/test_markov_regression.py index f98b9fbe0c6..0b187d12ae5 100644 --- a/statsmodels/tsa/regime_switching/tests/test_markov_regression.py +++ b/statsmodels/tsa/regime_switching/tests/test_markov_regression.py @@ -939,6 +939,86 @@ 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_py_hamilton_filter(self): mod = self.model params = self.true['params'] From d15b760628e0cc554e1bf48bd996e325671bf4f4 Mon Sep 17 00:00:00 2001 From: Gerard Brunick Date: Tue, 27 Feb 2018 20:43:31 -0500 Subject: [PATCH 2/3] Correct Hamilton filter docstrings * Give the correct shapes for all arguments. * Check that arguments shapes are consistent. * Check that transition matrices are left stochastic in the slow version of the function. This seems like a nice thing to do for the user, as most presentations use right stochastic matrices, so it is pretty easy to get this wrong. --- .../tsa/regime_switching/markov_switching.py | 44 +++++++++++++++++-- .../tests/test_markov_regression.py | 20 +++++++++ 2 files changed, 60 insertions(+), 4 deletions(-) diff --git a/statsmodels/tsa/regime_switching/markov_switching.py b/statsmodels/tsa/regime_switching/markov_switching.py index 9182faa4d40..a3b372de7b1 100644 --- a/statsmodels/tsa/regime_switching/markov_switching.py +++ b/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) @@ -221,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,). @@ -255,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 = ( diff --git a/statsmodels/tsa/regime_switching/tests/test_markov_regression.py b/statsmodels/tsa/regime_switching/tests/test_markov_regression.py index 0b187d12ae5..eb6d3b8e21b 100644 --- a/statsmodels/tsa/regime_switching/tests/test_markov_regression.py +++ b/statsmodels/tsa/regime_switching/tests/test_markov_regression.py @@ -1019,6 +1019,26 @@ def test_hamilton_filter_order_zero_with_tvtp(self): 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]: + try: + func(initial_probabilities, + regime_transition, + conditional_likelihoods) + except ValueError: + pass + else: + raise AssertionError('Did not raise ValueError.') + def test_py_hamilton_filter(self): mod = self.model params = self.true['params'] From 23b9bfb58b1a11b498391456227ebc63c62f8cde Mon Sep 17 00:00:00 2001 From: Gerard Brunick Date: Tue, 17 Apr 2018 20:42:04 -0400 Subject: [PATCH 3/3] Use assert_raises --- .../tsa/regime_switching/tests/test_markov_regression.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/statsmodels/tsa/regime_switching/tests/test_markov_regression.py b/statsmodels/tsa/regime_switching/tests/test_markov_regression.py index eb6d3b8e21b..3d47383bdd9 100644 --- a/statsmodels/tsa/regime_switching/tests/test_markov_regression.py +++ b/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__)) @@ -1030,14 +1030,10 @@ def test_hamilton_filter_shape_checks(self): for func in [markov_switching.py_hamilton_filter, markov_switching.cy_hamilton_filter]: - try: + with assert_raises(ValueError): func(initial_probabilities, regime_transition, conditional_likelihoods) - except ValueError: - pass - else: - raise AssertionError('Did not raise ValueError.') def test_py_hamilton_filter(self): mod = self.model