code: def _sim_prediction_bayes(self, h, simulations):
               
  
          sim_vector = np.zeros([simulations,h])
  
          for n in range(0,simulations):
  
              t_z = self.draw_latent_variables(nsims=1).T[0]
              sigma2, Y, scores = self._model(t_z)
              t_z = np.array([self.latent_variables.z_list[k].prior.transform(t_z[k]) for k in range(t_z.shape[0])])
  
              # Create arrays to iteratre over        
              sigma2_exp = sigma2.copy()
              scores_exp = scores.copy()
  
              # Loop over h time periods          
              for t in range(0,h):
                  new_value = t_z[0]
  
                  if self.q != 0:
                      for j in range(1,self.q+1):
                          new_value += t_z[j]*scores_exp[-j]
  
                  if self.p != 0:
                      for k in range(1,self.p+1):
                          new_value += t_z[k+self.q]*sigma2_exp[-k]  
  
                  sigma2_exp = np.append(sigma2_exp,[new_value]) # For indexing consistency
                  scores_exp = np.append(scores_exp,scores[np.random.randint(scores.shape[0])]) # expectation of score is zero
  
              sim_vector[n] = sigma2_exp[-h:]
  
          return np.transpose(sim_vector)
docs: Simulates a h-step ahead mean prediction
  
          Parameters
          ----------
          h : int
              How many steps ahead for the prediction
  
          simulations : int
              How many simulations to perform
  
          Returns
          ----------
          Matrix of simulations