Creates the structure of the model (model matrices etc)

        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)

def _model(self, beta):
        """ Creates the structure of the model (model matrices etc)

        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)
        """     

        Y = np.array(self.data[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])])
        bias = z[self.latent_variables.z_indices['Bias']['start']:self.latent_variables.z_indices['Bias']['end']+1]
        bias = np.reshape(bias, (-1, self.units))
        output_bias = z[self.latent_variables.z_indices['Output bias']['start']]
        input_weights = z[self.latent_variables.z_indices['Input weight']['start']:self.latent_variables.z_indices['Input weight']['end']+1]
        input_weights = np.reshape(input_weights, (-1, self.units))
        output_weights = z[self.latent_variables.z_indices['Output weight']['start']:self.latent_variables.z_indices['Output weight']['end']+1]

        # Construct neural network
        h = self.activation(self.X.T.dot(input_weights) + bias[0])

        if self.layers > 1:
            hidden_weights = z[self.latent_variables.z_indices['Hidden weight']['start']:self.latent_variables.z_indices['Hidden weight']['end']+1]
            hidden_weights = np.reshape(hidden_weights, (self.layers-1, self.units, -1))

            for k in range(0, self.layers-1):
                h = self.activation(h.dot(hidden_weights[k]) + bias[1+k])

        return h.dot(output_weights) + output_bias, Y