Simulates a h-step ahead mean prediction

        Parameters
        ----------
        sigma2 : np.array
            The past predicted values

        Y : np.array
            The past data

        scores : np.array
            The past scores

        h : int
            How many steps ahead for the prediction

        t_params : np.array
            A vector of (transformed) latent variables

        simulations : int
            How many simulations to perform

        Returns
        ----------
        Matrix of simulations

def _sim_prediction(self, sigma2, Y, scores, h, t_params, simulations):
        """ Simulates a h-step ahead mean prediction

        Parameters
        ----------
        sigma2 : np.array
            The past predicted values

        Y : np.array
            The past data

        scores : np.array
            The past scores

        h : int
            How many steps ahead for the prediction

        t_params : np.array
            A vector of (transformed) latent variables

        simulations : int
            How many simulations to perform

        Returns
        ----------
        Matrix of simulations
        """     

        sim_vector = np.zeros([simulations,h])

        for n in range(0,simulations):
            # 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_params[0]

                if self.q != 0:
                    for j in range(1,self.q+1):
                        new_value += t_params[j]*scores_exp[-j]

                if self.p != 0:
                    for k in range(1,self.p+1):
                        new_value += t_params[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)