code: def _mean_prediction(self, mu, Y, h, t_z):
               
  
          # Create arrays to iteratre over
          Y_exp = Y.copy()
          mu_exp = mu.copy()
  
          # Loop over h time periods          
          for t in range(0,h):
              new_value = t_z[0]
  
              if self.ar != 0:
                  for j in range(1, self.ar+1):
                      new_value += t_z[j]*Y_exp[-j]
  
              if self.ma != 0:
                  for k in range(1, self.ma+1):
                      if (k-1) >= t:
                          new_value += t_z[k+self.ar]*(Y_exp[-k]-self.link(mu_exp[-k]))
  
              if self.model_name2 == "Exponential":
                  Y_exp = np.append(Y_exp, [1.0/self.link(new_value)])
              else:
                  Y_exp = np.append(Y_exp, [self.link(new_value)])
  
              mu_exp = np.append(mu_exp,[0]) # For indexing consistency
          del mu_exp
          return Y_exp
docs: Creates a h-step ahead mean prediction
  
          This function is used for predict(). We have to iterate over the number
          of timepoints (h) that the user wants to predict, using as inputs the ARIMA
          parameters, past datapoints, and past predicted datapoints.
  
          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
  
          Returns
          ----------
          h-length np.array of mean predictions