code: def _construct_predict(self, beta, h):    
                       
  
          # 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)
docs: 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