Simulates a h-step ahead mean prediction

        Parameters
        ----------
        lmda : np.array
            The combined volatility component

        lmda_c : np.array
            The two volatility components

        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, lmda, lmda_c, Y, scores, h, t_params, simulations):
        """ Simulates a h-step ahead mean prediction

        Parameters
        ----------
        lmda : np.array
            The combined volatility component

        lmda_c : np.array
            The two volatility components

        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
            lmda_exp = lmda.copy()
            lmda_c_exp = [lmda_c.T[0],lmda_c.T[1]] 
            scores_exp = scores.copy()
            Y_exp = Y.copy()

            # Loop over h time periods          
            for t in range(0,h):
                new_value_comp = np.zeros(2)
                new_value = 0

                for comp in range(2):

                    if self.p != 0:
                        for j in range(self.p):
                            new_value_comp[comp] += t_params[1+j+(comp*(self.q+self.p))]*lmda_c_exp[comp][-j-1]

                    if self.q != 0:
                        for k in range(self.q):
                            new_value_comp[comp] += t_params[1+k+self.p+(comp*(self.q+self.p))]*scores_exp[-k-1]

                if self.leverage is True:
                    new_value_comp[1] += t_params[-3]*np.sign(-(Y_exp[-1]-t_params[-1]))*(scores_exp[t-1]+1)

                lmda_exp = np.append(lmda_exp,[new_value_comp.sum()+t_params[0]]) # For indexing consistency
                lmda_c_exp[0] = np.append(lmda_c_exp[0],new_value_comp[0])
                lmda_c_exp[1] = np.append(lmda_c_exp[1],new_value_comp[1])
                rnd_no = np.random.randint(scores.shape[0])
                scores_exp = np.append(scores_exp,scores[rnd_no]) # expectation of score is zero
                Y_exp = np.append(Y_exp,Y[rnd_no])
            sim_vector[n] = lmda_exp[-h:]

        return np.transpose(sim_vector)