code: def _compute_pred_intervals(model, forecasts, h, level):
      sigma = model["sigma2"]
      season_length = model["m"]
      pf = forecasts["mean"]
  
      model_type = model["components"]
      steps = steps = np.arange(1, h + 1)
      hm = np.floor((h - 1) / season_length)
      last_state = model["states"][-1]
  
      # error, trend, and seasonality type
      error = model_type[0]
      trend = model_type[1]
      seasonality = model_type[2]
      damped = model_type[3]
  
      # parameters
      alpha = model["par"][0]
      beta = model["par"][1]
      gamma = model["par"][2]
      phi = model["par"][3]
  
      exp1 = (
          alpha**2
          + alpha * beta * steps
          + (1 / 6) * beta**2 * steps * (2 * steps - 1)
      )
      exp2 = (beta * phi * steps) / (1 - phi) ** 2
      exp3 = 2 * alpha * (1 - phi) + beta * phi
      exp4 = (beta * phi * (1 - phi**steps)) / (
          (1 - phi) ** 2 * (1 - phi**2)
      )
      exp5 = 2 * alpha * (1 - phi**2) + beta * phi * (
          1 + 2 * phi - phi**steps
      )
  
      compute_intervals = True
      # Class 1 models
      if (
          error == "A"
          and trend == "N"
          and seasonality == "N"
          and damped == "N"
      ):
          # Model ANN
          sigmah = 1 + alpha**2 * (steps - 1)
          sigmah = sigma * sigmah
  
      elif (
          error == "A"
          and trend == "A"
          and seasonality == "N"
          and damped == "N"
      ):
          # Model AAN
          sigmah = 1 + (steps - 1) * exp1
          sigmah = sigma * sigmah
  
      elif (
          error == "A"
          and trend == "A"
          and seasonality == "N"
          and damped == "D"
      ):
          # Model AAdN
          sigmah = (
              1 + alpha**2 * (steps - 1) + exp2 * exp3 - exp4 * exp5
          )
          sigmah = sigma * sigmah
  
      elif (
          error == "A"
          and trend == "N"
          and seasonality == "A"
          and damped == "N"
      ):
          # Model ANA
          sigmah = (
              1
              + alpha**2 * (steps - 1)
              + gamma * hm * (2 * alpha + gamma)
          )
          sigmah = sigma * sigmah
  
      elif (
          error == "A"
          and trend == "A"
          and seasonality == "A"
          and damped == "N"
      ):
          # Model AAA
          exp6 = 2 * alpha + gamma + beta * season_length * (hm + 1)
          sigmah = 1 + (steps - 1) * exp1 + gamma * hm * exp6
          sigmah = sigma * sigmah
  
      elif (
          error == "A"
          and trend == "A"
          and seasonality == "A"
          and damped == "D"
      ):
          # Model AAdA
          exp7 = (2 * beta * gamma * phi) / (
              (1 - phi) * (1 - phi**season_length)
          )
          exp8 = hm * (
              1 - phi**season_length
          ) - phi**season_length * (1 - phi ** (season_length * hm))
          sigmah = (
              1
              + alpha**2 * (steps - 1)
              + exp2 * exp3
              - exp4 * exp5
              + gamma * hm * (2 * alpha + gamma)
              + exp7 * exp8
          )
          sigmah = sigma * sigmah
  
      # Class 2 models
      elif (
          error == "M"
          and trend == "N"
          and seasonality == "N"
          and damped == "N"
      ):
          # Model MNN
          cvals = np.full(h, alpha)
          sigmah = _compute_sigmah(pf, h, sigma, cvals)
  
      elif (
          error == "M"
          and trend == "A"
          and seasonality == "N"
          and damped == "N"
      ):
          # Model MAN
          cvals = alpha + beta * steps
          sigmah = _compute_sigmah(pf, h, sigma, cvals)
  
      elif (
          error == "M"
          and trend == "A"
          and seasonality == "N"
          and damped == "D"
      ):
          # Model MAdN
          cvals = np.full(h, np.nan)
          for k in range(1, h + 1):
              sum_phi = 0
              for j in range(1, k + 1):
                  sum_phi = sum_phi + phi**j
              cvals[k - 1] = alpha + beta * sum_phi
          sigmah = _compute_sigmah(pf, h, sigma, cvals)
  
      elif (
          error == "M"
          and trend == "N"
          and seasonality == "A"
          and damped == "N"
      ):
          # Model MNA
          dvals = np.zeros(h)
          for k in range(1, h + 1):
              val = k % season_length
              if val == 0:
                  dvals[k - 1] = 1
          cvals = alpha + gamma * dvals
          sigmah = _compute_sigmah(pf, h, sigma, cvals)
  
      elif (
          error == "M"
          and trend == "A"
          and seasonality == "A"
          and damped == "N"
      ):
          # Model MAA
          dvals = np.zeros(h)
          for k in range(1, h + 1):
              val = k % season_length
              if val == 0:
                  dvals[k - 1] = 1
          cvals = alpha * beta * steps + gamma * dvals
          sigmah = _compute_sigmah(pf, h, sigma, cvals)
  
      elif (
          error == "M"
          and trend == "A"
          and seasonality == "A"
          and damped == "D"
      ):
          # Model MAdA
          dvals = np.zeros(h)
          for k in range(1, h + 1):
              val = k % season_length
              if val == 0:
                  dvals[k - 1] = 1
          cvals = np.full(h, np.nan)
          for k in range(1, h + 1):
              sum_phi = 0
              for j in range(1, k + 1):
                  sum_phi = sum_phi + phi**j
              cvals[k - 1] = (
                  alpha + beta * sum_phi + gamma * dvals[k - 1]
              )
          sigmah = _compute_sigmah(pf, h, sigma, cvals)
  
      elif error == "M" and seasonality == "M":
          # Class 3 models
          sigmah = _class3models(
              h,
              sigma,
              last_state,
              season_length,
              error,
              trend,
              seasonality,
              damped,
              alpha,
              beta,
              gamma,
              phi,
          )
  
      else:
          # Classes 4 and 5 models
          np.random.seed(1)
          compute_intervals = False
          nsim = 5000
          y_path = np.zeros([nsim, h])
  
          if math.isnan(beta):
              beta = 0
          if math.isnan(gamma):
              gamma = 0
          if math.isnan(phi):
              phi = 0
  
          for k in range(nsim):
              e = np.random.normal(0, np.sqrt(sigma), h)
              yhat = np.zeros(h)
              etssimulate(
                  last_state,
                  season_length,
                  switch(error),
                  switch(trend),
                  switch(seasonality),
                  alpha,
                  beta,
                  gamma,
                  phi,
                  h,
                  yhat,
                  e,
              )
              y_path[k,] = yhat
  
          lower = np.quantile(
              y_path, 0.5 - np.array(level) / 200, axis=0
          )
          upper = np.quantile(
              y_path, 0.5 + np.array(level) / 200, axis=0
          )
          pi = {
              **{f"lo-{lv}": lower[i] for i, lv in enumerate(level)},
              **{f"hi-{lv}": upper[i] for i, lv in enumerate(level)},
          }
  
      if compute_intervals:
          pi = _calculate_intervals(
              forecasts, level=level, h=h, sigmah=np.sqrt(sigmah)
          )
  
      return pi
  
imports: import math
  from collections import namedtuple
  from typing import Tuple
  import numpy
  from numba import njit
  from numba.typed import List
  from statsmodels.tsa.seasonal import seasonal_decompose
  from .utils import _calculate_intervals, CACHE, NOGIL
label_desc: Use image node_id 25 for calling a global function with example usage: _compute_pred_intervals(model, forecasts, h, level) and returns: pi