Simulates a h-step ahead mean prediction

        Parameters
        ----------
        h : int
            How many steps ahead for the prediction

        simulations : int
            How many simulations to perform

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

def _sim_prediction_bayes(self, h, simulations):
        """ Simulates a h-step ahead mean prediction

        Parameters
        ----------
        h : int
            How many steps ahead for the prediction

        simulations : int
            How many simulations to perform

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

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

        for n in range(0,simulations):

            t_z = self.draw_latent_variables(nsims=1).T[0]
            theta, theta_t, Y, scores = self._model(t_z)
            t_z = np.array([self.latent_variables.z_list[k].prior.transform(t_z[k]) for k in range(t_z.shape[0])])

            model_scale, model_shape, model_skewness = self._get_scale_and_shape(t_z)

            Y_exp = Y.copy()
            theta_exp = theta.copy()
            theta_t_exp = theta_t.copy()
            scores_exp = scores.copy()

            #(TODO: vectorize the inner construction here)      
            for t in range(0,h):
                new_value1 = theta_t_exp[-1] + theta_exp[-1] + t_z[0]*scores_exp[-1]
                new_value2 = theta_t_exp[-1] + t_z[1]*scores_exp[-1]

                if self.model_name2 == "Exponential":
                    rnd_value = self.family.draw_variable(1.0/self.link(new_value1),model_scale,model_shape,model_skewness,1)[0]
                else:
                    rnd_value = self.family.draw_variable(self.link(new_value1),model_scale,model_shape,model_skewness,1)[0]

                Y_exp = np.append(Y_exp,[rnd_value])
                theta_exp = np.append(theta_exp,[new_value1]) # For indexing consistency
                theta_t_exp = np.append(theta_t_exp,[new_value2])
                scores_exp = np.append(scores_exp,scores[np.random.randint(scores.shape[0])]) # expectation of score is zero

            sim_vector[n] = Y_exp[-h:]

        return np.transpose(sim_vector)