Merge pull request #5383 from ChadFulton/exact-diffuse-sim

BUG: Initialization of simulation smoother with exact diffuse initialization
statsmodels · Nov 15, 2018 · 5c71c8e · 5c71c8e
2 parents 14fb572 + 24687d7
commit 5c71c8e
Show file tree

Hide file tree

Showing 2 changed files with 31 additions and 0 deletions.
diff --git a/statsmodels/tsa/statespace/_simulation_smoother.pyx.in b/statsmodels/tsa/statespace/_simulation_smoother.pyx.in
@@ -378,10 +378,12 @@ cdef class {{prefix}}SimulationSmoother(object):
             raise RuntimeError("Statespace model not initialized.")
         blas.{{prefix}}copy(&k_states, &self.model.initial_state[0], &inc, &self.simulated_model.initial_state[0], &inc)
         blas.{{prefix}}copy(&k_states2, &self.model.initial_state_cov[0,0], &inc, &self.simulated_model.initial_state_cov[0,0], &inc)
+        blas.{{prefix}}copy(&k_states2, &self.model.initial_diffuse_state_cov[0,0], &inc, &self.simulated_model.initial_diffuse_state_cov[0,0], &inc)
 
         if self.has_missing:
             blas.{{prefix}}copy(&k_states, &self.model.initial_state[0], &inc, &self.secondary_simulated_model.initial_state[0], &inc)
             blas.{{prefix}}copy(&k_states2, &self.model.initial_state_cov[0,0], &inc, &self.secondary_simulated_model.initial_state_cov[0,0], &inc)
+            blas.{{prefix}}copy(&k_states2, &self.model.initial_diffuse_state_cov[0,0], &inc, &self.secondary_simulated_model.initial_diffuse_state_cov[0,0], &inc)
 
         # 0. Kalman filter initialization: get alpha_1^+ ~ N(a_1, P_1)
         # Usually, this means transforming the N(0,1) random variate

diff --git a/statsmodels/tsa/statespace/tests/test_simulation_smoothing.py b/statsmodels/tsa/statespace/tests/test_simulation_smoothing.py
@@ -637,3 +637,32 @@ def test_simulation_smoothing_state_intercept():
     sim.simulate(disturbance_variates=np.zeros(mod.nobs * 2),
                  initial_state_variates=np.zeros(1))
     assert_equal(sim.simulated_state[0], intercept)
+
+
+def test_simulation_smoothing_state_intercept_diffuse():
+    nobs = 10
+    intercept = 100
+    endog = np.ones(nobs) * intercept
+
+    # Test without missing values
+    mod = sarimax.SARIMAX(endog, order=(0, 0, 0), trend='c',
+                          measurement_error=True,
+                          initialization='diffuse')
+    mod.update([intercept, 1., 1.])
+
+    sim = mod.simulation_smoother()
+    sim.simulate(disturbance_variates=np.zeros(mod.nobs * 2),
+                 initial_state_variates=np.zeros(1))
+    assert_equal(sim.simulated_state[0], intercept)
+
+    # Test with missing values
+    endog[5] = np.nan
+    mod = sarimax.SARIMAX(endog, order=(0, 0, 0), trend='c',
+                          measurement_error=True,
+                          initialization='diffuse')
+    mod.update([intercept, 1., 1.])
+
+    sim = mod.simulation_smoother()
+    sim.simulate(disturbance_variates=np.zeros(mod.nobs * 2),
+                 initial_state_variates=np.zeros(1))
+    assert_equal(sim.simulated_state[0], intercept)