code: def tbats_model(
      y,
      seasonal_periods,
      use_boxcox,
      bc_lower_bound,
      bc_upper_bound,
      use_trend,
      use_damped_trend,
      use_arma_errors,
  ):
      # Sort seasonal periods
      seasonal_periods = np.sort(seasonal_periods)
  
      # Check if there are missing values
      indices = np.where(np.isnan(y))[0]
      if len(indices) > 0:
          max_index = indices[-1]
          y = y[max_index + 1 : len(y)]
  
      # Check if there are negative values
      if np.any(y < 0):
          use_boxcox = False
  
      # Check if there is a trend
      if use_trend:
          adf = adfuller(y, regression="ct")
          if adf[1] <= 0.05:
              warnings.warn(
                  "The time series is trend-stationary, disabling trend components."
              )
              use_trend = False
              use_damped_trend = False
  
      # Choose the number of harmonics
      k_vector = np.array(
          [find_harmonics(y, m) for m in seasonal_periods]
      )
  
      # First model - No ARMA errors
      ar_coeffs = None
      ma_coeffs = None
  
      best_model = tbats_model_generator(
          y,
          seasonal_periods,
          k_vector,
          use_boxcox,
          bc_lower_bound,
          bc_upper_bound,
          use_trend,
          use_damped_trend,
          use_arma_errors,
          ar_coeffs,
          ma_coeffs,
      )
      best_model["description"] = "Model without ARMA errors"
  
      # Second model
      # Uses ARMA errors and the k_vector from the first model
      if use_arma_errors:
          errors = best_model["errors"][
              0
          ]  # ARMA errors from first model
          with np.errstate(invalid="ignore"):
              fit = auto_arima_f(errors, max_d=0, seasonal=False)
          p, q = fit["arma"][:2]
          if p != 0:
              ar_coeffs = np.zeros(p)
          else:
              ar_coeffs = None
          if q != 0:
              ma_coeffs = np.zeros(q)
          else:
              ma_coeffs = None
  
          model_arma_errors = tbats_model_generator(
              y,
              seasonal_periods,
              k_vector,
              use_boxcox,
              bc_lower_bound,
              bc_upper_bound,
              use_trend,
              use_damped_trend,
              use_arma_errors,
              ar_coeffs,
              ma_coeffs,
          )
  
          if model_arma_errors["aic"] <= best_model["aic"]:
              best_model = model_arma_errors
              best_model["description"] = "Model with ARMA errors"
  
      return best_model
  
imports: import warnings
  from functools import partial
  import numpy
  import pandas
  import scipy.linalg
  import statsmodels.api
  from numba import njit
  from numpy.polynomial.polynomial import Polynomial
  from scipy.special import inv_boxcox
  from scipy.stats import boxcox
  from scipy.optimize import minimize
  from statsmodels.tsa.stattools import adfuller
  from .arima import auto_arima_f
  from .utils import NOGIL, CACHE
label_desc: Use image node_id 16 for calling a global function with example usage: tbats_model(y, seasonal_periods, use_boxcox, bc_lower_bound, bc_upper_bound, use_trend, use_damped_trend, use_arma_errors) and returns: best_model