code: def _model(self, beta):
          
  
          Y = np.array(self.data[self.max_lag:self.data.shape[0]])
          X = np.ones(Y.shape[0])
          scores = np.zeros(Y.shape[0])
  
          # Transform latent variables
          parm = np.array([self.latent_variables.z_list[k].prior.transform(beta[k]) for k in range(beta.shape[0])])
          lmda = np.ones(Y.shape[0])*parm[0]
  
          # Loop over time series
          for t in range(0, Y.shape[0]):
  
              if t < self.max_lag:
  
                  lmda[t] = parm[0]/(1-np.sum(parm[1:(self.p+1)]))
  
              else:
  
                  # Loop over GARCH terms
                  for p_term in range(0, self.p):
                      lmda[t] += parm[1+p_term]*lmda[t-p_term-1]
  
                  # Loop over Score terms
                  for q_term in range(0, self.q):
                      lmda[t] += parm[1+self.p+q_term]*scores[t-q_term-1]
  
                  if self.leverage is True:
                      lmda[t] += parm[-3]*np.sign(-(Y[t-1]-parm[-1]))*(scores[t-1]+1)
  
              scores[t] = (((parm[-2]+1.0)*np.power(Y[t]-parm[-1],2))/float(parm[-2]*np.exp(lmda[t]) + np.power(Y[t]-parm[-1],2))) - 1.0
          return lmda, Y, scores
docs: Creates the structure of the model (model matrices, etc)
          
          Parameters
          ----------
          
          beta : np.array
              Contains untransformed starting values for the latent variables
          
          Returns
          ----------
          
          lambda : np.array
              Contains the values for the conditional volatility series
          
          Y : np.array
              Contains the length-adjusted time series (accounting for lags)
          
          scores : np.array
              Contains the score terms for the time series