Creates h-step ahead forecasts for the Gaussian process
        
        Parameters
        ----------
        beta : np.array
            Contains untransformed starting values for the latent variables
        
        h: int
            How many steps ahead to forecast
        
        Returns
        ----------
        - predictions
        - variance of predictions

def _construct_predict(self, beta, h):    
        """ Creates h-step ahead forecasts for the Gaussian process
        
        Parameters
        ----------
        beta : np.array
            Contains untransformed starting values for the latent variables
        
        h: int
            How many steps ahead to forecast
        
        Returns
        ----------
        - predictions
        - variance of predictions
        """             

        # Refactor this entire code in future
        parm = np.array([self.latent_variables.z_list[k].prior.transform(beta[k]) for k in range(beta.shape[0])])
        Xstart = self.X().copy()
        Xstart = [i for i in Xstart]
        predictions = np.zeros(h)
        variances = np.zeros(h)

        for step in range(0,h):
            Xstar = []

            for lag in range(0,self.max_lag):
                if lag == 0:
                    if step == 0:
                        Xstar.append([self.data[-1]])
                        Xstart[0] = np.append(Xstart[0],self.data[-1])
                    else:
                        Xstar.append([predictions[step-1]])
                        Xstart[0] = np.append(Xstart[0],predictions[step-1])
                else:
                    Xstar.append([Xstart[lag-1][-2]])
                    Xstart[lag] = np.append(Xstart[lag],Xstart[lag-1][-2])

            Kstar = self.kernel.Kstar(parm, np.transpose(np.array(Xstar)))

            L = self._L(parm)
            alpha = self._alpha(L)   

            predictions[step] = np.dot(np.transpose(Kstar), alpha)
            v = la.cho_solve((L, True), Kstar)
            variances[step] = self.kernel.Kstarstar(parm, np.transpose(np.array(Xstar))) - np.dot(v.T, v)

        return predictions, variances, predictions - 1.98*np.power(variances,0.5), predictions + 1.98*np.power(variances,0.5)