code: def _sim_prediction(self, mu, Y, h, t_z, simulations):
               
  
          model_scale, model_shape, model_skewness = self._get_scale_and_shape(t_z)
          sim_vector = np.zeros([simulations,h])
  
          for n in range(0, simulations):
              # Create arrays to iteratre over        
              Y_exp = Y.copy()
  
              # Loop over h time periods          
              for t in range(0,h):
  
                  if self.ar != 0:
                      Y_exp_normalized = (Y_exp[-self.ar:][::-1] - self._norm_mean) / self._norm_std
                      new_value = self.predict_new(np.append(1.0, Y_exp_normalized), self.latent_variables.get_z_values())
  
                  else:
                      new_value = self.predict_new(np.array([1.0]), self.latent_variables.get_z_values())
  
                  new_value += np.random.randn(1)*t_z[-1]
  
                  if self.model_name2 == "Exponential":
                      rnd_value = self.family.draw_variable(1.0/self.link(new_value), model_scale, model_shape, model_skewness, 1)[0]
                  else:
                      rnd_value = self.family.draw_variable(self.link(new_value), model_scale, model_shape, model_skewness, 1)[0]
  
                  Y_exp = np.append(Y_exp, [rnd_value])
  
              sim_vector[n] = Y_exp[-h:]
  
          return np.transpose(sim_vector)
docs: Simulates a h-step ahead mean prediction
  
          Parameters
          ----------
          mu : np.ndarray
              The past predicted values
  
          Y : np.ndarray
              The past data
  
          h : int
              How many steps ahead for the prediction
  
          t_z : np.ndarray
              A vector of (transformed) latent variables
  
          simulations : int
              How many simulations to perform
  
          Returns
          ----------
          Matrix of simulations