BUG: simulation smoothed measurement disturbance with FILTER_COLLAPSED

statsmodels/tsa/statespace/_simulation_smoother.pyx.in

@@ -32,8 +32,9 @@ cdef int SIMULATE_ALL = (
 )
 
 from statsmodels.tsa.statespace._kalman_filter cimport (
-    FILTER_CONVENTIONAL, INVERT_UNIVARIATE, SOLVE_CHOLESKY,
-    TIMING_INIT_PREDICTED, STABILITY_FORCE_SYMMETRY, MEMORY_STORE_ALL
+    FILTER_CONVENTIONAL, FILTER_COLLAPSED, FILTER_UNIVARIATE,
+    INVERT_UNIVARIATE, SOLVE_CHOLESKY, TIMING_INIT_PREDICTED,
+    STABILITY_FORCE_SYMMETRY, MEMORY_STORE_ALL
 )
 from statsmodels.tsa.statespace._kalman_smoother cimport (
     SMOOTHER_ALL
@@ -237,7 +238,10 @@ cdef class {{prefix}}SimulationSmoother(object):
 
         # Simulated data (\tilde eta_t, \tilde eps_t, \tilde alpha_t)
         # Note that these are (k_endog x nobs), (k_posdef x nobs), (k_states x nobs)
-        dim2[0] = self.model.k_endog; dim2[1] = self.nobs;
+        if filter_method & FILTER_COLLAPSED:
+            dim2[0] = self.model.k_states; dim2[1] = self.nobs;
+        else:
+            dim2[0] = self.model.k_endog; dim2[1] = self.nobs;
         self.simulated_measurement_disturbance = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN)
         dim2[0] = self.model.k_posdef; dim2[1] = self.nobs;
         self.simulated_state_disturbance = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN)
@@ -264,9 +268,9 @@ cdef class {{prefix}}SimulationSmoother(object):
         self._tmp2 = &self.tmp2[0,0]
 
     def __reduce__(self):
-        args = (self.model, self.filter_method, self.inversion_method,
-                self.stability_method, self.conserve_memory, self.filter_timing,
-                self.tolerance, self.loglikelihood_burn, self.smoother_output,
+        args = (self.model, self.model.filter_method, self.model.inversion_method,
+                self.model.stability_method, self.model.conserve_memory, self.model.filter_timing,
+                self.model.tolerance, self.model.loglikelihood_burn, self.simulated_smoother.smoother_output,
                 self.simulation_output, self.nobs, self.pretransformed_variates)
         state = {
             'disturbance_variates': np.array(self.disturbance_variates, copy=True, order='F'),
@@ -360,6 +364,9 @@ cdef class {{prefix}}SimulationSmoother(object):
         cdef:
             {{cython_type}} alpha = 1.0
             {{cython_type}} gamma = -1.0
+        cdef {{prefix}}Statespace collapsing_model
+        cdef {{cython_type}} [:,:] collapsing_obs
+        cdef np.npy_intp dim2[2]
 
 
         if simulation_output == -1:
@@ -485,6 +492,11 @@ cdef class {{prefix}}SimulationSmoother(object):
         self.simulated_smoother.smoother_output = simulation_output
         self.simulated_smoother()
 
+        # In the FILTER_COLLAPSED case, we now need to replace nobs_endog with
+        # nobs_kstates
+        if self.simulated_kfilter.filter_method & FILTER_COLLAPSED:
+            nobs_endog = nobs_kstates
+
         if self.has_missing:
             # This gives us \hat w_t
             #               \hat alpha_t+1
@@ -520,12 +532,34 @@ cdef class {{prefix}}SimulationSmoother(object):
         #               \tilde alpha_t+1 = \hat alpha_t^* + alpha_t^+  (simulation_output & SIMULATE_STATE)
         if self.simulation_output & SIMULATE_DISTURBANCE:
             # \tilde eps_t = \hat eps_t^* + eps_t^+
-            blas.{{prefix}}copy(&nobs_endog, &self.disturbance_variates[0], &inc, &self.simulated_measurement_disturbance[0,0], &inc)
+            # In the FILTER_COLLAPSED case, we have to use a collapsed version of the disturbance variates,
+            # so we construct a new representation that can perform the collapse on the variates
+            if self.simulated_kfilter.filter_method & FILTER_COLLAPSED:
+                dim2[0] = k_endog; dim2[1] = self.nobs;
+                collapsing_obs = np.PyArray_ZEROS(2, dim2, {{typenum}}, FORTRAN)
+                nobs_endog = k_endog * self.nobs
+                blas.{{prefix}}copy(&nobs_endog, &self.disturbance_variates[0], &inc, &collapsing_obs[0, 0], &inc)
+                nobs_endog = nobs_kstates
+
+                collapsing_model = {{prefix}}Statespace(
+                    collapsing_obs,
+                    self.model.design, self.model.obs_intercept, self.model.obs_cov,
+                    self.model.transition, self.model.state_intercept, self.model.selection,
+                    self.model.state_cov
+                )
+                collapsing_model.initialized = True
+                for t in range(self.nobs):
+                    collapsing_model.seek(t, self.simulated_kfilter.filter_method & FILTER_UNIVARIATE,
+                                             self.simulated_kfilter.filter_method & FILTER_COLLAPSED)
+                    blas.{{prefix}}copy(&k_states, &collapsing_model.collapse_obs[0], &inc, &self.simulated_measurement_disturbance[0,t], &inc)
+            else:
+                blas.{{prefix}}copy(&nobs_endog, &self.disturbance_variates[0], &inc, &self.simulated_measurement_disturbance[0,0], &inc)
+
             blas.{{prefix}}axpy(&nobs_endog, &alpha, &self.simulated_smoother.smoothed_measurement_disturbance[0,0], &inc,
                                                      &self.simulated_measurement_disturbance[0,0], &inc)
 
             # \tilde eta_t = \hat eta_t^* + eta_t^+
-            blas.{{prefix}}copy(&nobs_posdef, &self.disturbance_variates[nobs_endog], &inc, &self.simulated_state_disturbance[0,0], &inc)
+            blas.{{prefix}}copy(&nobs_posdef, &self.disturbance_variates[k_endog * self.nobs], &inc, &self.simulated_state_disturbance[0,0], &inc)
             blas.{{prefix}}axpy(&nobs_posdef, &alpha, &self.simulated_smoother.smoothed_state_disturbance[0,0], &inc,
                                                      &self.simulated_state_disturbance[0,0], &inc)
 

statsmodels/tsa/statespace/simulation_smoother.py

@@ -522,7 +522,7 @@ def simulated_state_disturbance(self):
         return self._simulated_state_disturbance
 
     def simulate(self, simulation_output=-1, disturbance_variates=None,
-                 initial_state_variates=None):
+                 initial_state_variates=None, pretransformed_variates=False):
         r"""
         Perform simulation smoothing
 
@@ -564,7 +564,8 @@ def simulate(self, simulation_output=-1, disturbance_variates=None,
         # simulation
         if disturbance_variates is not None:
             self._simulation_smoother.set_disturbance_variates(
-                np.array(disturbance_variates, dtype=self.dtype)
+                np.array(disturbance_variates, dtype=self.dtype),
+                pretransformed=pretransformed_variates
             )
         else:
             self._simulation_smoother.draw_disturbance_variates()
@@ -573,7 +574,8 @@ def simulate(self, simulation_output=-1, disturbance_variates=None,
         # simulation
         if initial_state_variates is not None:
             self._simulation_smoother.set_initial_state_variates(
-                np.array(initial_state_variates, dtype=self.dtype)
+                np.array(initial_state_variates, dtype=self.dtype),
+                pretransformed=pretransformed_variates
             )
         else:
             self._simulation_smoother.draw_initial_state_variates()

statsmodels/tsa/statespace/tests/test_collapsed.py

@@ -418,6 +418,35 @@ def test_using_alterate(self):
 
 
 class TestDFM(object):
+    @classmethod
+    def create_model(cls, obs, **kwargs):
+        # Create the model with typical state space
+        mod = MLEModel(obs, k_states=2, k_posdef=2, **kwargs)
+        mod['design'] = np.array([[-32.47143586, 17.33779024],
+                                  [-7.40264169, 1.69279859],
+                                  [-209.04702853, 125.2879374]])
+        mod['obs_cov'] = np.diag(
+            np.array([0.0622668, 1.95666886, 58.37473642]))
+        mod['transition'] = np.array([[0.29935707, 0.33289005],
+                                      [-0.7639868, 1.2844237]])
+        mod['selection'] = np.eye(2)
+        mod['state_cov'] = np.array([[1.2, -0.25],
+                                     [-0.25, 1.1]])
+        mod.initialize_approximate_diffuse(1e6)
+        return mod
+
+    @classmethod
+    def collapse(cls, obs, **kwargs):
+        mod = cls.create_model(obs, **kwargs)
+        mod.smooth([], return_ssm=True)
+
+        _ss = mod.ssm._statespace
+        out = np.zeros((mod.nobs, mod.k_states))
+        for t in range(mod.nobs):
+            _ss.seek(t, mod.ssm.filter_univariate, 1)
+            out[t] = np.array(_ss.collapse_obs)
+        return out
+
     @classmethod
     def setup_class(cls, which='mixed', *args, **kwargs):
         # Data
@@ -438,40 +467,27 @@ def setup_class(cls, which='mixed', *args, **kwargs):
             obs.iloc[119:130, 0] = np.nan
             obs.iloc[119:130, 2] = np.nan
 
-        # Create the model with typical state space
-        mod = MLEModel(obs, k_states=2, k_posdef=2, **kwargs)
-        mod['design'] = np.array([[-32.47143586, 17.33779024],
-                                  [-7.40264169, 1.69279859],
-                                  [-209.04702853, 125.2879374]])
-        mod['obs_cov'] = np.diag(
-            np.array([0.0622668, 1.95666886, 58.37473642]))
-        mod['transition'] = np.array([[0.29935707, 0.33289005],
-                                      [-0.7639868, 1.2844237]])
-        mod['selection'] = np.eye(2)
-        mod['state_cov'] = np.array([[1.2, -0.25],
-                                     [-0.25, 1.1]])
-        mod.initialize_approximate_diffuse(1e6)
+        mod = cls.create_model(obs, **kwargs)
         cls.model = mod.ssm
 
         n_disturbance_variates = (
             (cls.model.k_endog + cls.model.k_posdef) * cls.model.nobs
         )
+        np.random.seed(1234)
+        dv = np.random.normal(size=n_disturbance_variates)
+        isv = np.random.normal(size=cls.model.k_states)
 
         # Collapsed filtering, smoothing, and simulation smoothing
         cls.model.filter_collapsed = True
         cls.results_b = cls.model.smooth()
-        cls.sim_b = cls.model.simulation_smoother(
-            disturbance_variates=np.zeros(n_disturbance_variates),
-            initial_state_variates=np.zeros(cls.model.k_states)
-        )
+        cls.sim_b = cls.model.simulation_smoother()
+        cls.sim_b.simulate(disturbance_variates=dv, initial_state_variates=isv)
 
         # Conventional filtering, smoothing, and simulation smoothing
         cls.model.filter_collapsed = False
         cls.results_a = cls.model.smooth()
-        cls.sim_a = cls.model.simulation_smoother(
-            disturbance_variates=np.zeros(n_disturbance_variates),
-            initial_state_variates=np.zeros(cls.model.k_states)
-        )
+        cls.sim_a = cls.model.simulation_smoother()
+        cls.sim_a.simulate(disturbance_variates=dv, initial_state_variates=isv)
 
         # Create the model with augmented state space
         kwargs.pop('filter_collapsed', None)
@@ -618,19 +634,22 @@ def test_smoothed_state_disturbance_cov(self):
     def test_simulation_smoothed_state(self):
         assert_allclose(
             self.sim_a.simulated_state,
-            self.sim_a.simulated_state
+            self.sim_b.simulated_state
         )
 
+    # Skipped because "measurement" refers to different things; even different
+    # dimensions
+    @skip
     def test_simulation_smoothed_measurement_disturbance(self):
         assert_allclose(
             self.sim_a.simulated_measurement_disturbance,
-            self.sim_a.simulated_measurement_disturbance
+            self.sim_b.simulated_measurement_disturbance
         )
 
     def test_simulation_smoothed_state_disturbance(self):
         assert_allclose(
             self.sim_a.simulated_state_disturbance,
-            self.sim_a.simulated_state_disturbance
+            self.sim_b.simulated_state_disturbance
         )
 
 
@@ -675,15 +694,28 @@ def test_smooth_method(self):
                      SMOOTH_ALTERNATIVE)
 
 
-class TestDFMClassicalSmoothing(TestDFM):
+class TestDFMMeasurementDisturbance(TestDFM):
     @classmethod
     def setup_class(cls, *args, **kwargs):
-        super(TestDFMClassicalSmoothing, cls).setup_class(
-            smooth_method=SMOOTH_CLASSICAL, **kwargs)
+        super(TestDFMMeasurementDisturbance, cls).setup_class(
+            smooth_method=SMOOTH_CLASSICAL, which='none', **kwargs)
 
-    def test_smooth_method(self):
-        assert_equal(self.model.smooth_method, SMOOTH_CLASSICAL)
-        assert_equal(self.model._kalman_smoother.smooth_method,
-                     SMOOTH_CLASSICAL)
-        assert_equal(self.model._kalman_smoother._smooth_method,
-                     SMOOTH_CLASSICAL)
+    def test_smoothed_state_disturbance(self):
+        assert_allclose(
+            self.results_a.smoothed_state_disturbance,
+            self.results_b.smoothed_state_disturbance, atol=1e-7)
+
+    def test_smoothed_measurement_disturbance(self):
+        assert_allclose(
+            self.collapse(self.results_a.smoothed_measurement_disturbance.T).T,
+            self.results_b.smoothed_measurement_disturbance, atol=1e-7)
+
+    def test_simulation_smoothed_measurement_disturbance(self):
+        assert_allclose(
+            self.collapse(self.sim_a.simulated_measurement_disturbance.T),
+            self.sim_b.simulated_measurement_disturbance.T, atol=1e-7)
+
+    def test_simulation_smoothed_state_disturbance(self):
+        assert_allclose(
+            self.sim_a.simulated_state_disturbance,
+            self.sim_b.simulated_state_disturbance, atol=1e-7)