code: def forecast_theta(obj, h, level=None):
      forecast = np.full(h, fill_value=np.nan)
      n = obj["n"]
      states = obj["states"]
      alpha = obj["par"]["alpha"]
      theta = obj["par"]["theta"]
      thetaforecast(
          states=states,
          n=n,
          modeltype=switch_theta(obj["modeltype"]),
          h=h,
          f=forecast,
          alpha=alpha,
          theta=theta,
      )
      res = {"mean": forecast}
  
      if level is not None:
          sigma = np.std(obj["residuals"][3:], ddof=1)
          mean_y = obj["mean_y"]
          samples = compute_pi_samples(
              n=n,
              h=h,
              states=states,
              sigma=sigma,
              alpha=alpha,
              theta=theta,
              mean_y=mean_y,
          )
          for lv in level:
              min_q = (100 - lv) / 200
              max_q = min_q + lv / 100
              res[f"lo-{lv}"] = np.quantile(samples, min_q, axis=1)
              res[f"hi-{lv}"] = np.quantile(samples, max_q, axis=1)
  
      if obj.get("decompose", False):
          seas_forecast = _repeat_val_seas(
              obj["seas_forecast"]["mean"], h=h
          )
          for key in res:
              if obj["decomposition_type"] == "multiplicative":
                  res[key] = res[key] * seas_forecast
              else:
                  res[key] = res[key] + seas_forecast
      return res
  
imports: import math
  from typing import Tuple
  import numpy
  from numba import njit
  from scipy.stats import norm
  from statsmodels.tsa.seasonal import seasonal_decompose
  from statsmodels.tsa.stattools import acf
  from .ets import restrict_to_bounds, results
  from .utils import _seasonal_naive, _repeat_val_seas, CACHE, NOGIL
label_desc: Use image node_id 15 for calling a global function with example usage: forecast_theta(obj, h, level) and returns: res