code: def auto_arima_f(
      x,
      d=None,
      D=None,
      max_p=5,
      max_q=5,
      max_P=2,
      max_Q=2,
      max_order=5,
      max_d=2,
      max_D=1,
      start_p=2,
      start_q=2,
      start_P=1,
      start_Q=1,
      stationary=False,
      seasonal=True,
      ic="aicc",
      stepwise=True,
      nmodels=94,
      trace=False,
      approximation=None,
      method=None,
      truncate=None,
      xreg=None,
      test="kpss",
      test_kwargs=None,
      seasonal_test="seas",
      seasonal_test_kwargs=None,
      allowdrift=True,
      allowmean=True,
      blambda=None,
      biasadj=False,
      period=1,
  ):
      if approximation is None:
          approximation = len(x) > 150 or period > 12
      if x.ndim > 1:
          raise ValueError(
              "auto_arima can only handle univariate time series"
          )
      if test_kwargs is None:
          test_kwargs = {}
      if seasonal_test_kwargs is None:
          seasonal_test_kwargs = {}
      x = x.copy()
      origx = x
      missing = np.isnan(x)
      nonmissing_idxs = np.where(~missing)[0]
      firstnonmiss = nonmissing_idxs.min()
      lastnonmiss = nonmissing_idxs.max() + 1
      series_len = int(np.sum(~missing[firstnonmiss:lastnonmiss]))
      x = x[firstnonmiss:]
      if xreg is not None:
          if xreg.dtype not in (np.float32, np.float64):
              raise ValueError("xreg should be a float array")
          xreg = xreg[firstnonmiss:]
      if is_constant(x):
          if np.isnan(x).all():
              raise ValueError("all data are missing")
          if allowmean:
              fit = Arima(
                  x, order=(0, 0, 0), fixed=np.array([np.mean(x)])
              )
          else:
              fit = Arima(x, order=(0, 0, 0), include_mean=False)
          fit["x"] = origx
          fit["constant"] = True
          return fit
      m = period if seasonal else 1
      if m < 1:
          m = 1
      else:
          m = round(m)
      max_p = min(max_p, series_len // 3)
      max_q = min(max_q, series_len // 3)
      max_P = min(max_P, math.floor(series_len / 3 / m))
      max_Q = min(max_Q, math.floor(series_len / 3 / m))
      if series_len <= 3:
          ic = "aic"
      if blambda is not None:
          raise NotImplementedError("blambda != None")
          # x = boxcox(x, blambda)
          # setattr(blambda, 'biasadj', biasadj)
      if xreg is not None:
          xx = x.copy()
          xregg = xreg.copy()
          constant_columns = np.array(
              [is_constant(col) for col in xregg.T]
          )
          if constant_columns.all():
              xregg = None
          else:
              if constant_columns.any():
                  xregg = xregg[:, ~constant_columns]
              X = np.hstack(
                  [
                      np.arange(1, xregg.shape[0] + 1).reshape(-1, 1),
                      xregg,
                  ]
              )
              X = X[~np.isnan(X).any(1)]
              _, sv, _ = np.linalg.svd(X)
              if sv.min() / sv.sum() < np.finfo(np.float64).eps:
                  raise ValueError("xreg is rank deficient")
              j = (~np.isnan(x)) & (~np.isnan(np.nansum(xregg, 1)))
              xx[j] = sm.OLS(x, sm.add_constant(xregg)).fit().resid
      else:
          xx = x
          xregg = None
      if stationary:
          d = D = 0
      if m == 1:
          D = max_P = max_Q = 0
      elif D is None and len(xx) <= 2 * m:
          D = 0
      elif D is None:
          D = nsdiffs(
              xx,
              period=m,
              test=seasonal_test,
              max_D=max_D,
              **seasonal_test_kwargs,
          )
          if D > 0 and xregg is not None:
              diffxreg = diff(xregg, m, D)
              if any(is_constant(col) for col in xregg.T):
                  D -= 1
          if D > 0:
              dx = diff(xx, m, D)
              if np.isnan(dx).all():
                  D -= 1
      if D > 0:
          dx = diff(xx, m, D)
      else:
          dx = xx
      if xregg is not None:
          if D > 0:
              diffxreg = diff(xregg, m, D)
          else:
              diffxreg = xregg
      if d is None:
          d = ndiffs(dx, test=test, max_d=max_d, **test_kwargs)
          if d > 0 and xregg is not None:
              diffxreg = diff(diffxreg, 1, d)
              if any(is_constant(col) for col in diffxreg.T):
                  d -= 1
          if d > 0:
              diffdx = diff(dx, 1, d)
              if np.isnan(diffdx).all():
                  d -= 1
      if D >= 2:
          warnings.warn(
              "Having more than one seasonal differences is not recommended. Please consider using only one seasonal difference."
          )
      elif D + d > 2:
          warnings.warn(
              "Having 3 or more differencing operations is not recommended. Please consider reducing the total number of differences."
          )
      if d > 0:
          dx = diff(dx, 1, d)
      if len(dx) == 0:
          raise ValueError("not enough data to proceed")
      elif is_constant(dx):
          if xreg is None:
              if D > 0 and d == 0:
                  fit = Arima(
                      x,
                      order=(0, d, 0),
                      seasonal={"order": (0, D, 0), "period": m},
                      include_constant=True,
                      fixed=np.array([np.mean(dx / m)]),
                      method=method,
                  )
              elif D > 0 and d > 0:
                  fit = Arima(
                      x,
                      order=(0, d, 0),
                      seasonal={"order": (0, D, 0), "period": m},
                      method=method,
                  )
              elif d == 2:
                  fit = Arima(x, order=(0, d, 0), method=method)
              elif d < 2:
                  fit = Arima(
                      x,
                      order=(0, d, 0),
                      include_constant=True,
                      fixed=np.array([np.mean(dx)]),
                      method=method,
                  )
              else:
                  raise ValueError(
                      "Data follow a simple polynomial and are not suitable for ARIMA modelling."
                  )
          else:
              if D > 0:
                  fit = Arima(
                      x,
                      order=(0, d, 0),
                      seasonal={"order": (0, D, 0), "period": m},
                      xreg=xreg,
                      method=method,
                  )
              else:
                  fit = Arima(
                      x, order=(0, d, 0), xreg=xreg, method=method
                  )
          fit["x"] = origx
          return fit
      if m > 1:
          if max_p > 0:
              max_p = min(max_p, m - 1)
          if max_q > 0:
              max_q = min(max_q, m - 1)
      if approximation:
          if truncate is not None:
              if len(x) > truncate:
                  x = x[-truncate:]
          try:
              if D == 0:
                  fit = arima(x, order=(0, d, 0), xreg=xreg)
              else:
                  fit = arima(
                      x,
                      order=(0, d, 0),
                      seasonal={"order": (0, D, 0), "period": m},
                      xreg=xreg,
                  )
              offset = -2 * fit["loglik"] - series_len * math.log(
                  fit["sigma2"]
              )
          except:
              offset = 0
      else:
          offset = 0
      allowdrift = allowdrift and (d + D) == 1
      allowmean = allowmean and (d + D) == 0
      constant = allowdrift or allowmean
      if approximation and trace:
          print(
              "Fitting models using approximations to speed things up"
          )
      if not stepwise:
          bestfit = search_arima(
              x,
              d,
              D,
              max_p,
              max_q,
              max_P,
              max_Q,
              max_order,
              stationary,
              ic,
              trace,
              approximation,
              method=method,
              xreg=xreg,
              offset=offset,
              allowdrift=allowdrift,
              allowmean=allowmean,
              period=m,
          )
          bestfit["lambda"] = blambda
          bestfit["x"] = origx
          if trace:
              print(
                  f"Best model: {arima_string(bestfit, padding=True)}\n\n"
              )
          return bestfit
      if len(x) < 10:
          start_p = min(start_p, 1)
          start_q = min(start_q, 1)
          start_P = 0
          start_Q = 0
      p = start_p = min(start_p, max_p)
      q = start_q = min(start_q, max_q)
      P = start_P = min(start_P, max_P)
      Q = start_Q = min(start_Q, max_Q)
      results = np.full((nmodels, 8), np.nan)
      p_myarima = partial(
          myarima,
          x=x,
          constant=constant,
          ic=ic,
          trace=trace,
          approximation=approximation,
          offset=offset,
          xreg=xreg,
          method=method,
      )
      bestfit = p_myarima(
          order=(p, d, q),
          seasonal={"order": (P, D, Q), "period": m},
      )
      results[0] = (p, d, q, P, D, Q, constant, bestfit["ic"])
      fit = p_myarima(
          order=(0, d, 0),
          seasonal={"order": (0, D, 0), "period": m},
      )
      results[1] = (0, d, 0, 0, D, 0, constant, fit["ic"])
      if fit["ic"] < bestfit["ic"]:
          bestfit = fit
          p = q = P = Q = 0
      k = 1
      if max_p > 0 or max_P > 0:
          p_ = int(max_p > 0)
          P_ = int(m > 1 and max_P > 0)
          fit = p_myarima(
              order=(p_, d, 0),
              seasonal={"order": (P_, D, 0), "period": m},
          )
          results[k + 1] = (p_, d, 0, P_, D, 0, constant, fit["ic"])
          if fit["ic"] < bestfit["ic"]:
              bestfit = fit
              p = p_
              P = P_
              q = Q = 0
          k += 1
      if max_q > 0 or max_Q > 0:
          q_ = int(max_q > 0)
          Q_ = int(m > 1 and max_Q > 0)
          fit = p_myarima(
              order=(0, d, q_),
              seasonal={"order": (0, D, Q_), "period": m},
          )
          results[k + 1] = (0, d, q_, 0, D, Q_, constant, fit["ic"])
          if fit["ic"] < bestfit["ic"]:
              bestfit = fit
              p = P = 0
              Q = Q_
              q = q_
          k += 1
      if constant:
          fit = p_myarima(
              order=(0, d, 0),
              seasonal={"order": (0, D, 0), "period": m},
              constant=False,
          )
          results[k + 1] = (0, d, 0, 0, D, 0, 0, fit["ic"])
          if fit["ic"] < bestfit["ic"]:
              bestfit = fit
              p = q = P = Q = 0
          k += 1
  
      def try_params(p, d, q, P, D, Q, constant, k, bestfit):
          k += 1
          improved = False
          if k >= results.shape[0]:
              return k, bestfit, improved
          fit = p_myarima(
              order=(p, d, q),
              seasonal={"order": (P, D, Q), "period": m},
          )
          results[k] = (p, d, q, P, D, Q, constant, fit["ic"])
          if fit["ic"] < bestfit["ic"]:
              bestfit = fit
              improved = True
          return k, bestfit, improved
  
      startk = 0
      while startk < k and k < nmodels:
          startk = k
          if P > 0 and newmodel(
              p, d, q, P - 1, D, Q, constant, results[:k]
          ):
              k, bestfit, improved = try_params(
                  p, d, q, P - 1, D, Q, constant, k, bestfit
              )
              if improved:
                  P -= 1
                  continue
          if Q > 0 and newmodel(
              p, d, q, P, D, Q - 1, constant, results[:k]
          ):
              k, bestfit, improved = try_params(
                  p, d, q, P, D, Q - 1, constant, k, bestfit
              )
              if improved:
                  Q -= 1
                  continue
          if P < max_P and newmodel(
              p, d, q, P + 1, D, Q, constant, results[:k]
          ):
              k, bestfit, improved = try_params(
                  p, d, q, P + 1, D, Q, constant, k, bestfit
              )
              if improved:
                  P += 1
                  continue
          if Q < max_Q and newmodel(
              p, d, q, P, D, Q + 1, constant, results[:k]
          ):
              k, bestfit, improved = try_params(
                  p, d, q, P, D, Q + 1, constant, k, bestfit
              )
              if improved:
                  Q += 1
                  continue
          if (
              Q > 0
              and P > 0
              and newmodel(
                  p, d, q, P - 1, D, Q - 1, constant, results[:k]
              )
          ):
              k, bestfit, improved = try_params(
                  p, d, q, P - 1, D, Q - 1, constant, k, bestfit
              )
              if improved:
                  P -= 1
                  Q -= 1
                  continue
          if (
              Q < max_Q
              and P > 0
              and newmodel(
                  p, d, q, P - 1, D, Q + 1, constant, results[:k]
              )
          ):
              k, bestfit, improved = try_params(
                  p, d, q, P - 1, D, Q + 1, constant, k, bestfit
              )
              if improved:
                  P -= 1
                  Q += 1
                  continue
          if (
              Q > 0
              and P < max_P
              and newmodel(
                  p, d, q, P + 1, D, Q - 1, constant, results[:k]
              )
          ):
              k, bestfit, improved = try_params(
                  p, d, q, P + 1, D, Q - 1, constant, k, bestfit
              )
              if improved:
                  P += 1
                  Q -= 1
                  continue
          if (
              Q < max_Q
              and P < max_P
              and newmodel(
                  p, d, q, P + 1, D, Q + 1, constant, results[:k]
              )
          ):
              k, bestfit, improved = try_params(
                  p, d, q, P + 1, D, Q + 1, constant, k, bestfit
              )
              if improved:
                  P += 1
                  Q += 1
                  continue
          if p > 0 and newmodel(
              p - 1, d, q, P, D, Q, constant, results[:k]
          ):
              k, bestfit, improved = try_params(
                  p - 1, d, q, P, D, Q, constant, k, bestfit
              )
              if improved:
                  p -= 1
                  continue
          if q > 0 and newmodel(
              p, d, q - 1, P, D, Q, constant, results[:k]
          ):
              k, bestfit, improved = try_params(
                  p, d, q - 1, P, D, Q, constant, k, bestfit
              )
              if improved:
                  q -= 1
                  continue
          if p < max_p and newmodel(
              p + 1, d, q, P, D, Q, constant, results[:k]
          ):
              k, bestfit, improved = try_params(
                  p + 1, d, q, P, D, Q, constant, k, bestfit
              )
              if improved:
                  p += 1
                  continue
          if q < max_q and newmodel(
              p, d, q + 1, P, D, Q, constant, results[:k]
          ):
              k, bestfit, improved = try_params(
                  p, d, q + 1, P, D, Q, constant, k, bestfit
              )
              if improved:
                  q += 1
                  continue
          if (
              q > 0
              and p > 0
              and newmodel(
                  p - 1, d, q - 1, P, D, Q, constant, results[:k]
              )
          ):
              k, bestfit, improved = try_params(
                  p - 1, d, q - 1, P, D, Q, constant, k, bestfit
              )
              if improved:
                  p -= 1
                  q -= 1
                  continue
          if (
              q < max_q
              and p > 0
              and newmodel(
                  p - 1, d, q + 1, P, D, Q, constant, results[:k]
              )
          ):
              k, bestfit, improved = try_params(
                  p - 1, d, q + 1, P, D, Q, constant, k, bestfit
              )
              if improved:
                  p -= 1
                  q += 1
                  continue
          if (
              q > 0
              and p < max_p
              and newmodel(
                  p + 1, d, q - 1, P, D, Q, constant, results[:k]
              )
          ):
              k, bestfit, improved = try_params(
                  p + 1, d, q - 1, P, D, Q, constant, k, bestfit
              )
              if improved:
                  p += 1
                  q -= 1
                  continue
          if (
              q < max_q
              and p < max_p
              and newmodel(
                  p + 1, d, q + 1, P, D, Q, constant, results[:k]
              )
          ):
              k, bestfit, improved = try_params(
                  p + 1, d, q + 1, P, D, Q, constant, k, bestfit
              )
              if improved:
                  p += 1
                  q += 1
                  continue
          if (allowdrift or allowmean) and newmodel(
              p, d, q, P, D, Q, not constant, results[:k]
          ):
              k, bestfit, improved = try_params(
                  p, d, q, P, D, Q, not constant, k, bestfit
              )
              if improved:
                  constant = not constant
                  continue
      if k >= nmodels:
          warnings.warn(
              f"Stepwise search was stopped early due to reaching the model number limit: nmodels={nmodels}"
          )
      if approximation or bestfit["arma"] is not None:
          if trace:
              print(
                  "Now re-fitting the best model(s) without approximations...\n"
              )
          icorder = np.argsort(results[:, 7])
          nmodels = np.sum(~np.isnan(results[:, 7]))
          for i in range(nmodels):
              k = icorder[i]
              p, q, P, Q, constant = map(
                  int, results[k, [0, 2, 3, 5, 6]]
              )
              fit = myarima(
                  x,
                  (p, d, q),
                  {"order": (P, D, Q), "period": m},
                  constant=results[k, 6],
                  ic=ic,
                  trace=trace,
                  approximation=False,
                  method=method,
                  xreg=xreg,
              )
              if fit["ic"] < math.inf:
                  bestfit = fit
                  break
      if math.isinf(bestfit["ic"]) and method != "CSS":
          raise ValueError("No suitable ARIMA model found")
  
      bestfit["x"] = origx
      bestfit["ic"] = None
      bestfit["lambda"] = blambda
  
      return bestfit
  
imports: import math
  import warnings
  from collections import namedtuple
  from functools import partial
  from typing import Optional, Dict, Union, Tuple
  import numpy
  import pandas
  import statsmodels.api
  from numba import njit
  from scipy.optimize import minimize
  from scipy.stats import norm
  from .mstl import mstl
  from .utils import CACHE, NOGIL
label_desc: Use image node_id 36 for calling a global function with example usage: auto_arima_f(x, d, D, max_p, max_q, max_P, max_Q, max_order, max_d, max_D, start_p, start_q, start_P, start_Q, stationary, seasonal, ic, stepwise, nmodels, trace, approximation, method, truncate, xreg, test, test_kwargs, seasonal_test, seasonal_test_kwargs, allowdrift, allowmean, blambda, biasadj, period) and returns: bestfit, fit, bestfit, k,  bestfit,  improved, fit, k,  bestfit,  improved