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

def _mean_prediction(self, mu, Y, h, t_z):
        """ 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
        """     

        # 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