Repository: statsmodels/statsmodels (Python)

Issue:
STLForecastResults.get_prediction() with ETSModel as model fails when "trend"="mul" in the model kwargs
I get an error when calling the get_prediction() method on a STLForecastResults object that is the output of STLForecast(endog=data, model=ETSModel, model_kwargs={'trend': 'mul'}).fit()

https://github.com/statsmodels/statsmodels/blob/ecd6ecaf09e4864af156bd549002700e84097377/statsmodels/tsa/exponential_smoothing/ets.py#L2290
`self.forecast_variance` is not defined when `self.method == "simulated"`
but `get_prediction()` calls `pred.var_pred_mean` that uses `self.forecast_variance`
https://github.com/statsmodels/statsmodels/blob/ecd6ecaf09e4864af156bd549002700e84097377/statsmodels/tsa/forecasting/stl.py#L504


Hints:
Stlf_model = STLForecast( ts, model =ARIMA, model_kwargs ={ 'order': (2,1,0)}, period =52)

Stlf_model = Stlf_model.fit()

Stlf_model.get_prediction()   gives the following error 

' Data must be 1-dimensional'

However, Stlf_model.forecast() runs fine.

But I need to get the corresponding confidence interval of the forecasted values
Can you provide a complete, copy-pastable example that reproduces?
Hi, 

I am using the following code (data attached)


import pandas as pd 
from statsmodels.tsa.forecasting.stl import STLForecast
from statsmodels.tsa.arima.model import ARIMA


ts  = pd. read_csv('../ts_example.csv')
ts  = ts.set_index('Date')
ts  = ts.resample('W-SUN').interpolate(method='linear')
ts  = ts.sort_index()  


stlf_model    = STLForecast(endog=ts,model = ARIMA,model_kwargs= { 'order':(2,1,0)}, period = 52)
stlf_model    = stlf_model.fit()
forecast      = stlf_model.forecast(26)

stlf_model.get_prediction() —- throws an error



[ts_example.zip](https://github.com/statsmodels/statsmodels/files/6076016/ts_example.zip)





Patch:
diff --git a/statsmodels/tsa/exponential_smoothing/base.py b/statsmodels/tsa/exponential_smoothing/base.py
--- a/statsmodels/tsa/exponential_smoothing/base.py
+++ b/statsmodels/tsa/exponential_smoothing/base.py
@@ -206,7 +206,10 @@ def _wrap_data(self, data, start_idx, end_idx, names=None):
         if data.ndim > 1 and data.shape[1] == 1:
             data = np.squeeze(data, axis=1)
         if self.use_pandas:
-            _, _, _, index = self._get_prediction_index(start_idx, end_idx)
+            if data.shape[0]:
+                _, _, _, index = self._get_prediction_index(start_idx, end_idx)
+            else:
+                index = None
             if data.ndim < 2:
                 data = pd.Series(data, index=index, name=names)
             else:
diff --git a/statsmodels/tsa/exponential_smoothing/ets.py b/statsmodels/tsa/exponential_smoothing/ets.py
--- a/statsmodels/tsa/exponential_smoothing/ets.py
+++ b/statsmodels/tsa/exponential_smoothing/ets.py
@@ -1757,7 +1757,6 @@ def simulate(
            principles and practice*, 2nd edition, OTexts: Melbourne,
            Australia. OTexts.com/fpp2. Accessed on February 28th 2020.
         """
-
         # Get the starting location
         start_idx = self._get_prediction_start_index(anchor)
 
@@ -2286,20 +2285,20 @@ def __init__(
                     )
                     # anchor
                     anchor = start_smooth + i
-            sim_results.append(
-                results.simulate(
-                    ndynamic,
-                    anchor=anchor_dynamic,
-                    repetitions=simulate_repetitions,
-                    **simulate_kwargs,
+            if ndynamic:
+                sim_results.append(
+                    results.simulate(
+                        ndynamic,
+                        anchor=anchor_dynamic,
+                        repetitions=simulate_repetitions,
+                        **simulate_kwargs,
+                    )
                 )
-            )
-            self.simulation_results = np.concatenate(sim_results, axis=0)
-            # if self.use_pandas:
-            #     self.simulation_results = pd.DataFrame(
-            #         self.simulation_results, index=self.row_labels,
-            #         columns=sim_results[0].columns
-            #     )
+            if sim_results and isinstance(sim_results[0], pd.DataFrame):
+                self.simulation_results = pd.concat(sim_results, 0)
+            else:
+                self.simulation_results = np.concatenate(sim_results, axis=0)
+            self.forecast_variance = self.simulation_results.var(1)
         else:  # method == 'exact'
             steps = np.ones(ndynamic + nsmooth)
             if ndynamic > 0:
diff --git a/statsmodels/tsa/forecasting/stl.py b/statsmodels/tsa/forecasting/stl.py
--- a/statsmodels/tsa/forecasting/stl.py
+++ b/statsmodels/tsa/forecasting/stl.py
@@ -498,8 +498,22 @@ def get_prediction(
         pred = self._model_result.get_prediction(
             start=start, end=end, dynamic=dynamic, **kwargs
         )
-        seasonal_prediction = self._get_seasonal_prediction(start, end, dynamic)
+        seasonal_prediction = self._get_seasonal_prediction(
+            start, end, dynamic
+        )
         mean = pred.predicted_mean + seasonal_prediction
+        try:
+            var_pred_mean = pred.var_pred_mean
+        except (AttributeError, NotImplementedError):
+            # Allow models that do not return var_pred_mean
+            import warnings
+
+            warnings.warn(
+                f"The variance of the predicted mean is not available using "
+                f"the {self.model.__class__.__name__} model class.",
+                UserWarning,
+            )
+            var_pred_mean = np.nan + mean.copy()
         return PredictionResults(
-            mean, pred.var_pred_mean, dist="norm", row_labels=pred.row_labels
+            mean, var_pred_mean, dist="norm", row_labels=pred.row_labels
         )


Test patch:
diff --git a/statsmodels/tsa/forecasting/tests/test_stl.py b/statsmodels/tsa/forecasting/tests/test_stl.py
--- a/statsmodels/tsa/forecasting/tests/test_stl.py
+++ b/statsmodels/tsa/forecasting/tests/test_stl.py
@@ -3,7 +3,7 @@
 import pandas as pd
 import pytest
 
-from statsmodels.datasets import sunspots
+import statsmodels.datasets
 from statsmodels.tsa.ar_model import AutoReg
 from statsmodels.tsa.arima.model import ARIMA
 from statsmodels.tsa.base.prediction import PredictionResults
@@ -133,15 +133,58 @@ def fit(self, *args, **kwargs):
         STLForecast(data, FakeModelSummary).fit().summary()
 
 
-def test_get_prediction():
-    # GH7309
-    df = sunspots.load_pandas().data
+@pytest.fixture(scope="function")
+def sunspots():
+    df = statsmodels.datasets.sunspots.load_pandas().data
     df.index = np.arange(df.shape[0])
-    y = df.iloc[:, 0]
+    return df.iloc[:, 0]
+
+
+def test_get_prediction(sunspots):
+    # GH7309
     stlf_model = STLForecast(
-        y, model=ARIMA, model_kwargs={"order": (2, 2, 0)}, period=11
+        sunspots, model=ARIMA, model_kwargs={"order": (2, 2, 0)}, period=11
     )
     stlf_res = stlf_model.fit()
     pred = stlf_res.get_prediction()
     assert pred.predicted_mean.shape == (309,)
     assert pred.var_pred_mean.shape == (309,)
+
+
+@pytest.mark.parametrize("not_implemented", [True, False])
+def test_no_var_pred(sunspots, not_implemented):
+    class DummyPred:
+        def __init__(self, predicted_mean, row_labels):
+            self.predicted_mean = predicted_mean
+            self.row_labels = row_labels
+
+            def f():
+                raise NotImplementedError
+
+            if not_implemented:
+                self.forecast = property(f)
+
+    class DummyRes:
+        def __init__(self, res):
+            self._res = res
+
+        def forecast(self, *args, **kwargs):
+            return self._res.forecast(*args, **kwargs)
+
+        def get_prediction(self, *args, **kwargs):
+            pred = self._res.get_prediction(*args, **kwargs)
+
+            return DummyPred(pred.predicted_mean, pred.row_labels)
+
+    class DummyMod:
+        def __init__(self, y):
+            self._mod = ARIMA(y)
+
+        def fit(self, *args, **kwargs):
+            res = self._mod.fit(*args, **kwargs)
+            return DummyRes(res)
+
+    stl_mod = STLForecast(sunspots, model=DummyMod, period=11)
+    stl_res = stl_mod.fit()
+    pred = stl_res.get_prediction()
+    assert np.all(np.isnan(pred.var_pred_mean))