code: def simulation_smoother(self,beta):
                   
  
          T, Z, R, Q, H = self._ss_matrices(beta)
  
          # Generate e_t+ and n_t+
          rnd_h = np.random.normal(0,np.sqrt(H),self.data.shape[0]+1)
          q_dist = ss.multivariate_normal([0.0], Q)
          rnd_q = q_dist.rvs(self.data.shape[0]+1)
  
          # Generate a_t+ and y_t+
          a_plus = np.zeros((T.shape[0],self.data.shape[0]+1)) 
          a_plus[0,0] = np.mean(self.data[0:5])
          y_plus = np.zeros(self.data.shape[0])
  
          for t in range(0,self.data.shape[0]+1):
              if t == 0:
                  a_plus[:,t] = np.dot(T,a_plus[:,t]) + rnd_q[t]
                  y_plus[t] = np.dot(Z,a_plus[:,t]) + rnd_h[t]
              else:
                  if t != self.data.shape[0]:
                      a_plus[:,t] = np.dot(T,a_plus[:,t-1]) + rnd_q[t]
                      y_plus[t] = np.dot(Z,a_plus[:,t]) + rnd_h[t]
  
          alpha_hat,_ = self.smoothed_state(self.data,beta)
          alpha_hat_plus,_ = self.smoothed_state(y_plus,beta)
          alpha_tilde = alpha_hat - alpha_hat_plus + a_plus
      
          return alpha_tilde
docs: Koopman's simulation smoother - simulates from states given
          model latent variables and observations
  
          Parameters
          ----------
  
          beta : np.array
              Contains untransformed starting values for latent variables
  
          Returns
          ----------
          - A simulated state evolution