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

def _model(self, beta):
        """ 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
        """

        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