title: Simulating Time Series With Unobserved Components Model
body: <p>After fitting a local level model using <code>UnobservedComponents</code> from <code>statsmodels</code> , we are trying to find ways to simulate new time series with the results. Something like:</p>
  
  <pre><code>import numpy as np
  import statsmodels as sm
  from statsmodels.tsa.statespace.structural import UnobservedComponents
  
  np.random.seed(12345)
  ar = np.r_[1, 0.9]
  ma = np.array([1])
  arma_process = sm.tsa.arima_process.ArmaProcess(ar, ma)
  
  X = 100 + arma_process.generate_sample(nsample=100)
  y = 1.2 * x + np.random.normal(size=100)
  y[70:] += 10
  
  plt.plot(X, label='X')
  plt.plot(y, label='y')
  plt.axvline(69, linestyle='--', color='k')
  plt.legend();
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/soB8U.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/soB8U.png" alt="time series example"></a></p>
  
  <pre><code>ss = {}
  ss["endog"] = y[:70]
  ss["level"] = "llevel"
  ss["exog"] = X[:70]
  
  model = UnobservedComponents(**ss)
  trained_model = model.fit()
  </code></pre>
  
  <p>Is it possible to use <code>trained_model</code> to simulate new time series given the exogenous variable <code>X[70:]</code>? Just as we have the <code>arma_process.generate_sample(nsample=100)</code>, we were wondering if we could do something like:</p>
  
  <pre><code>trained_model.generate_random_series(nsample=100, exog=X[70:])
  </code></pre>
  
  <p>The motivation behind it is so that we can compute the probability of having a time series as extreme as the observed <code>y[70:]</code> (p-value for identifying the response is bigger than the predicted one).</p>
  
  <p><strong>[EDIT]</strong></p>
  
  <p>After reading Josef's and cfulton's comments, I tried implementing the following:</p>
  
  <pre><code>mod1 = UnobservedComponents(np.zeros(y_post), 'llevel', exog=X_post)
  mod1.simulate(f_model.params, len(X_post))
  </code></pre>
  
  <p>But this resulted in simulations that doesn't seem to track the <code>predicted_mean</code> of the forecast for <code>X_post</code> as exog. Here's an example:</p>
  
  <p><a href="https://i.stack.imgur.com/mLajH.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/mLajH.png" alt="enter image description here"></a></p>
  
  <p>While the <code>y_post</code> meanders around 100, the simulation is at -400. This approach always leads to p_value of 50%. </p>
  
  <p>So when I tried using the <code>initial_sate=0</code> and the random shocks, here's the result:</p>
  
  <p><a href="https://i.stack.imgur.com/hkhb0.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/hkhb0.png" alt="enter image description here"></a></p>
  
  <p>It seemed now that the simulations were following the predicted mean and its 95% credible interval (as cfulton commented below, this is actually a wrong approach as well as it's replacing the level variance of the trained model). </p>
  
  <p>I tried using this approach just to see what p-values I'd observe. Here's how I compute the p-value:</p>
  
  <pre><code>samples = 1000
  r = 0
  y_post_sum = y_post.sum()
  for _ in range(samples):
      sim = mod1.simulate(f_model.params, len(X_post), initial_state=0, state_shocks=np.random.normal(size=len(X_post)))
      r += sim.sum() &gt;= y_post_sum
  print(r / samples)
  </code></pre>
  
  <p>For context, this is the <a href="https://google.github.io/CausalImpact/CausalImpact.html" rel="nofollow noreferrer">Causal Impact</a> model developed by Google. As it's been implemented in R, we've been trying to replicate the implementation in Python using <code>statsmodels</code> as the core to process time series.</p>
  
  <p>We already have a quite cool WIP <a href="https://github.com/jamalsenouci/causalimpact" rel="nofollow noreferrer">implementation</a> but we still need to have the p-value to know when in fact we had an impact that is not explained by mere randomness (the approach of simulating series and counting the ones whose summation surpasses <code>y_post.sum()</code> is also implemented in Google's model).</p>
  
  <p>In my example I used <code>y[70:]</code> += 10. If I add just one instead of ten, Google's p-value computation returns <code>0.001</code> (there's an impact in <code>y</code>) whereas in Python's approach it's returning 0.247 (no impact).</p>
  
  <p>Only when I add +5 to <code>y_post</code> is that the model returns p_value of 0.02 and as it's lower than 0.05, we consider that there's an impact in <code>y_post</code>.</p>
  
  <p>I'm using python3, statsmodels version 0.9.0</p>
  
  <p><strong>[EDIT2]</strong></p>
  
  <p>After reading cfulton's comments I decided to fully debug the code to see what was happening. Here's what I found:</p>
  
  <p>When we create an object of type <code>UnobservedComponents</code>, eventually the representation of the Kalman Filter is initiated. As default, it <a href="https://github.com/statsmodels/statsmodels/blob/v0.9.0/statsmodels/tsa/statespace/representation.py#L105" rel="nofollow noreferrer">receives the parameter <code>initial_variance</code></a> as 1e6 which sets the <a href="https://github.com/statsmodels/statsmodels/blob/v0.9.0/statsmodels/tsa/statespace/representation.py#L327" rel="nofollow noreferrer">same property</a> of the object.</p>
  
  <p>When we run the <code>simulate</code> method, the <code>initial_state_cov</code> value <a href="https://github.com/statsmodels/statsmodels/blob/v0.9.0/statsmodels/tsa/statespace/structural.py#L779" rel="nofollow noreferrer">is created</a> using this same value:</p>
  
  <pre><code>initial_state_cov = (
          np.eye(self.k_states, dtype=self.ssm.transition.dtype) *
          self.ssm.initial_variance
      )
  </code></pre>
  
  <p>Finally, this same value is used to find <a href="https://github.com/statsmodels/statsmodels/blob/v0.9.0/statsmodels/tsa/statespace/kalman_filter.py#L953" rel="nofollow noreferrer">initial_state</a>:</p>
  
  <pre><code>initial_state = np.random.multivariate_normal(
      self._initial_state, self._initial_state_cov)
  </code></pre>
  
  <p>Which results in a normal distribution with 1e6 of standard deviation.</p>
  
  <p>I tried running the following then:</p>
  
  <pre><code>mod1 = UnobservedComponents(np.zeros(len(X_post)), level='llevel', exog=X_post, initial_variance=1)
  sim = mod1.simulate(f_model.params, len(X_post))
  plt.plot(sim, label='simul')
  plt.plot(y_post, label='y')
  plt.legend();
  print(sim.sum() &gt; y_post.sum())
  </code></pre>
  
  <p>Which resulted in:</p>
  
  <p><a href="https://i.stack.imgur.com/skfys.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/skfys.png" alt="enter image description here"></a></p>
  
  <p>I tested then the p-value and finally for a variation of +1 in <code>y_post</code> the model now is identifying correctly the added signal.</p>
  
  <p>Still, when I tested with the same data that we have in R's Google package the p-value was still off. Maybe it's a matter of further tweaking the input to increase its accuracy.</p>