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