code: def _normal_model(self, beta):
               
  
          Y = self.y[self.max_lag:]
  
          # Transform latent variables
          z = np.array([self.latent_variables.z_list[k].prior.transform(beta[k]) for k in range(beta.shape[0])])
  
          # Constant and AR terms
          if self.ar == 0:
              mu = np.transpose(self.ar_matrix)
          elif self.ar == 1:
              mu = np.transpose(self.ar_matrix)*z[:-self.family_z_no-self.ma-len(self.X_names)][0]
          else:
              mu = np.matmul(np.transpose(self.ar_matrix),z[:-self.family_z_no-self.ma-len(self.X_names)])
  
          # X terms
          mu = mu + np.matmul(self.X[self.integ+self.max_lag:],z[self.ma+self.ar:(self.ma+self.ar+len(self.X_names))])
  
          # MA terms
          if self.ma != 0:
              mu = arimax_recursion(z, mu, Y, self.max_lag, Y.shape[0], self.ar, self.ma)
  
          return mu, Y
docs: Creates the structure of the model (model matrices, etc) for
          a Normal family ARIMAX model.
  
          Parameters
          ----------
          beta : np.ndarray
              Contains untransformed starting values for the latent variables
  
          Returns
          ----------
          mu : np.ndarray
              Contains the predicted values (location) for the time series
  
          Y : np.ndarray
              Contains the length-adjusted time series (accounting for lags)