Repository: statsmodels/statsmodels

Issue:
initial_level doesn't appear to be respected in ExponentialSmoothing class
#### Describe the bug

[A clear and concise description of what the bug is. This should explain **why** the current behaviour is a problem and why the expected output is a better solution.]

#### Code Sample, a copy-pastable example if possible


```python
from statsmodels.tsa.holtwinters import ExponentialSmoothing
series = [0.0, 0.0, 0.0, 100.0, 0.0, 0.0, 0.0]
es = ExponentialSmoothing(
    series,
    initialization_method='known',
    initial_level=20.0
)
es_fit = es.fit()
es_fit.params
```

Outputs:
```python
{'smoothing_level': 1.4901161193847656e-08,
 'smoothing_trend': nan,
 'smoothing_seasonal': nan,
 'damping_trend': nan,
 'initial_level': 14.239682618110056,
 'initial_trend': nan,
 'initial_seasons': array([], dtype=float64),
 'use_boxcox': False,
 'lamda': None,
 'remove_bias': False}
```


<details>

**Note**: As you can see, there are many issues on our GitHub tracker, so it is very possible that your issue has been posted before. Please check first before submitting so that we do not have to handle and close duplicates.

**Note**: Please be sure you are using the latest released version of `statsmodels`, or a recent build of `main`. If your problem has been fixed in an unreleased version, you might be able to use `main` until a new release occurs. 

**Note**: If you are using a released version, have you verified that the bug exists in the main branch of this repository? It helps the limited resources if we know problems exist in the current main branch so that they do not need to check whether the code sample produces a bug in the next release.

</details>


If the issue has not been resolved, please file it in the issue tracker.

#### Expected Output

`es_fit.params` should return 20.0 in the `initial_level` field

Additionally, when calling `es_fit.fittedvalues` I was expecting the first returned value to be 20.0

#### Output of ``import statsmodels.api as sm; sm.show_versions()``

<details>

[paste the output of ``import statsmodels.api as sm; sm.show_versions()`` here below this line]

INSTALLED VERSIONS
------------------
Python: 3.10.6.final.0
OS: Linux 5.10.147-133.644.amzn2.x86_64 #1 SMP Fri Oct 14 01:16:24 UTC 2022 x86_64
byteorder: little
LC_ALL: None
LANG: en_US.UTF-8

statsmodels
===========

Installed: 0.13.2 (/home/ec2-user/anaconda3/envs/python3/lib/python3.10/site-packages/statsmodels)

Required Dependencies
=====================

cython: 0.29.32 (/home/ec2-user/anaconda3/envs/python3/lib/python3.10/site-packages/Cython)
numpy: 1.22.4 (/home/ec2-user/anaconda3/envs/python3/lib/python3.10/site-packages/numpy)
scipy: 1.9.3 (/home/ec2-user/anaconda3/envs/python3/lib/python3.10/site-packages/scipy)
pandas: 1.5.1 (/home/ec2-user/anaconda3/envs/python3/lib/python3.10/site-packages/pandas)
    dateutil: 2.8.2 (/home/ec2-user/anaconda3/envs/python3/lib/python3.10/site-packages/dateutil)
patsy: 0.5.3 (/home/ec2-user/anaconda3/envs/python3/lib/python3.10/site-packages/patsy)

Optional Dependencies
=====================

matplotlib: 3.6.1 (/home/ec2-user/anaconda3/envs/python3/lib/python3.10/site-packages/matplotlib)
    backend: module://matplotlib_inline.backend_inline 
cvxopt: Not installed
joblib: 1.2.0 (/home/ec2-user/anaconda3/envs/python3/lib/python3.10/site-packages/joblib)

Developer Tools
================

IPython: 7.33.0 (/home/ec2-user/anaconda3/envs/python3/lib/python3.10/site-packages/IPython)
    jinja2: 3.1.2 (/home/ec2-user/anaconda3/envs/python3/lib/python3.10/site-packages/jinja2)
sphinx: 5.3.0 (/home/ec2-user/anaconda3/envs/python3/lib/python3.10/site-packages/sphinx)
    pygments: 2.13.0 (/home/ec2-user/anaconda3/envs/python3/lib/python3.10/site-packages/pygments)
pytest: 7.2.0 (/home/ec2-user/anaconda3/envs/python3/lib/python3.10/site-packages/pytest)
virtualenv: Not installed

</details>


Discussion:

Patch:
diff --git a/statsmodels/tsa/holtwinters/_smoothers.py b/statsmodels/tsa/holtwinters/_smoothers.py
index 4b10a5998..3f198cd07 100644
--- a/statsmodels/tsa/holtwinters/_smoothers.py
+++ b/statsmodels/tsa/holtwinters/_smoothers.py
@@ -192,7 +192,7 @@ def holt_mul_dam(x, hw_args: HoltWintersArgs):
     for i in range(1, hw_args.n):
         lvl[i] = (y_alpha[i - 1]) + (alphac * (lvl[i - 1] * b[i - 1] ** phi))
         b[i] = (beta * (lvl[i] / lvl[i - 1])) + (betac * b[i - 1] ** phi)
-    return hw_args.y - lvl * b ** phi
+    return hw_args.y - lvl * b**phi
 
 
 def holt_add_dam(x, hw_args: HoltWintersArgs):
@@ -337,7 +337,7 @@ def holt_win_mul_mul_dam(x, hw_args: HoltWintersArgs):
         s[i + m - 1] = (y_gamma[i - 1] / (lvl[i - 1] * b[i - 1] ** phi)) + (
             gammac * s[i - 1]
         )
-    return hw_args.y - (lvl * b ** phi) * s[: -(m - 1)]
+    return hw_args.y - (lvl * b**phi) * s[: -(m - 1)]
 
 
 def holt_win_add_add_dam(x, hw_args: HoltWintersArgs):
diff --git a/statsmodels/tsa/holtwinters/model.py b/statsmodels/tsa/holtwinters/model.py
index 0f2eea308..1e8d885a6 100644
--- a/statsmodels/tsa/holtwinters/model.py
+++ b/statsmodels/tsa/holtwinters/model.py
@@ -287,8 +287,8 @@ class ExponentialSmoothing(TimeSeriesModel):
         )
         estimated = self._initialization_method == "estimated"
         self._estimate_level = estimated
-        self._estimate_trend = estimated and self.trend
-        self._estimate_seasonal = estimated and self.seasonal
+        self._estimate_trend = estimated and self.trend is not None
+        self._estimate_seasonal = estimated and self.seasonal is not None
         self._bounds = self._check_bounds(bounds)
         self._use_boxcox = use_boxcox
         self._lambda = np.nan
@@ -765,8 +765,6 @@ class ExponentialSmoothing(TimeSeriesModel):
         beta = data.beta
         phi = data.phi
         gamma = data.gamma
-        initial_level = data.level
-        initial_trend = data.trend
         y = data.y
         start_params = data.params
 
@@ -796,11 +794,11 @@ class ExponentialSmoothing(TimeSeriesModel):
                 alpha is None,
                 has_trend and beta is None,
                 has_seasonal and gamma is None,
-                initial_level is None,
-                has_trend and initial_trend is None,
+                self._estimate_level,
+                self._estimate_trend,
                 damped_trend and phi is None,
             ]
-            + [has_seasonal] * m,
+            + [has_seasonal and self._estimate_seasonal] * m,
         )
         (
             sel,
diff --git a/statsmodels/tsa/holtwinters/results.py b/statsmodels/tsa/holtwinters/results.py
index e5baa1dd6..07b9f5147 100644
--- a/statsmodels/tsa/holtwinters/results.py
+++ b/statsmodels/tsa/holtwinters/results.py
@@ -668,7 +668,7 @@ class HoltWintersResults(Results):
             resid = self.model._y - fitted
         else:
             resid = (self.model._y - fitted) / fitted
-        sigma = np.sqrt(np.sum(resid ** 2) / (len(resid) - n_params))
+        sigma = np.sqrt(np.sum(resid**2) / (len(resid) - n_params))
 
         # get random error eps
         if isinstance(random_errors, np.ndarray):


Test patch:
diff --git a/statsmodels/tsa/holtwinters/tests/test_holtwinters.py b/statsmodels/tsa/holtwinters/tests/test_holtwinters.py
index 17f134652..0e8f28b13 100644
--- a/statsmodels/tsa/holtwinters/tests/test_holtwinters.py
+++ b/statsmodels/tsa/holtwinters/tests/test_holtwinters.py
@@ -43,6 +43,62 @@ housing_data = housing_data.asfreq("MS")
 SEASONALS = ("add", "mul", None)
 TRENDS = ("add", "mul", None)
 
+# aust = pd.read_json(aust_json, typ='Series').sort_index()
+data = [
+    41.727457999999999,
+    24.04185,
+    32.328102999999999,
+    37.328707999999999,
+    46.213152999999998,
+    29.346326000000001,
+    36.482909999999997,
+    42.977719,
+    48.901524999999999,
+    31.180221,
+    37.717880999999998,
+    40.420211000000002,
+    51.206862999999998,
+    31.887228,
+    40.978262999999998,
+    43.772491000000002,
+    55.558566999999996,
+    33.850915000000001,
+    42.076383,
+    45.642291999999998,
+    59.766779999999997,
+    35.191876999999998,
+    44.319737000000003,
+    47.913736,
+]
+index = [
+    "2005-03-01 00:00:00",
+    "2005-06-01 00:00:00",
+    "2005-09-01 00:00:00",
+    "2005-12-01 00:00:00",
+    "2006-03-01 00:00:00",
+    "2006-06-01 00:00:00",
+    "2006-09-01 00:00:00",
+    "2006-12-01 00:00:00",
+    "2007-03-01 00:00:00",
+    "2007-06-01 00:00:00",
+    "2007-09-01 00:00:00",
+    "2007-12-01 00:00:00",
+    "2008-03-01 00:00:00",
+    "2008-06-01 00:00:00",
+    "2008-09-01 00:00:00",
+    "2008-12-01 00:00:00",
+    "2009-03-01 00:00:00",
+    "2009-06-01 00:00:00",
+    "2009-09-01 00:00:00",
+    "2009-12-01 00:00:00",
+    "2010-03-01 00:00:00",
+    "2010-06-01 00:00:00",
+    "2010-09-01 00:00:00",
+    "2010-12-01 00:00:00",
+]
+idx = pd.to_datetime(index)
+aust = pd.Series(data, index=pd.DatetimeIndex(idx, freq=pd.infer_freq(idx)))
+
 
 @pytest.fixture(scope="module")
 def ses():
@@ -240,63 +296,6 @@ class TestHoltWinters:
         )
         cls.livestock2_livestock = livestock2_livestock
 
-        # aust = pd.read_json(aust_json, typ='Series').sort_index()
-        data = [
-            41.727457999999999,
-            24.04185,
-            32.328102999999999,
-            37.328707999999999,
-            46.213152999999998,
-            29.346326000000001,
-            36.482909999999997,
-            42.977719,
-            48.901524999999999,
-            31.180221,
-            37.717880999999998,
-            40.420211000000002,
-            51.206862999999998,
-            31.887228,
-            40.978262999999998,
-            43.772491000000002,
-            55.558566999999996,
-            33.850915000000001,
-            42.076383,
-            45.642291999999998,
-            59.766779999999997,
-            35.191876999999998,
-            44.319737000000003,
-            47.913736,
-        ]
-        index = [
-            "2005-03-01 00:00:00",
-            "2005-06-01 00:00:00",
-            "2005-09-01 00:00:00",
-            "2005-12-01 00:00:00",
-            "2006-03-01 00:00:00",
-            "2006-06-01 00:00:00",
-            "2006-09-01 00:00:00",
-            "2006-12-01 00:00:00",
-            "2007-03-01 00:00:00",
-            "2007-06-01 00:00:00",
-            "2007-09-01 00:00:00",
-            "2007-12-01 00:00:00",
-            "2008-03-01 00:00:00",
-            "2008-06-01 00:00:00",
-            "2008-09-01 00:00:00",
-            "2008-12-01 00:00:00",
-            "2009-03-01 00:00:00",
-            "2009-06-01 00:00:00",
-            "2009-09-01 00:00:00",
-            "2009-12-01 00:00:00",
-            "2010-03-01 00:00:00",
-            "2010-06-01 00:00:00",
-            "2010-09-01 00:00:00",
-            "2010-12-01 00:00:00",
-        ]
-        aust = pd.Series(data, index)
-        aust.index = pd.DatetimeIndex(
-            aust.index, freq=pd.infer_freq(aust.index)
-        )
         cls.aust = aust
         cls.start_params = [
             1.5520372162082909e-09,
@@ -519,7 +518,11 @@ class TestHoltWinters:
         # livestock2_livestock <- c(...)
         # res <- ets(livestock2_livestock, model='AAN', damped_trend=TRUE,
         #            phi=0.98)
-        mod = Holt(self.livestock2_livestock, damped_trend=True)
+        mod = Holt(
+            self.livestock2_livestock,
+            damped_trend=True,
+            initialization_method="estimated",
+        )
         params = {
             "smoothing_level": 0.97402626,
             "smoothing_trend": 0.00010006,
@@ -1646,11 +1649,11 @@ def test_error_boxcox():
     with pytest.raises(TypeError, match="use_boxcox must be True"):
         ExponentialSmoothing(y, use_boxcox="a", initialization_method="known")
 
-    mod = ExponentialSmoothing(y ** 2, use_boxcox=True)
+    mod = ExponentialSmoothing(y**2, use_boxcox=True)
     assert isinstance(mod, ExponentialSmoothing)
 
     mod = ExponentialSmoothing(
-        y ** 2, use_boxcox=True, initialization_method="legacy-heuristic"
+        y**2, use_boxcox=True, initialization_method="legacy-heuristic"
     )
     with pytest.raises(ValueError, match="use_boxcox was set"):
         mod.fit(use_boxcox=False)
@@ -1950,7 +1953,7 @@ def test_attributes(ses):
 
 def test_summary_boxcox(ses):
     mod = ExponentialSmoothing(
-        ses ** 2, use_boxcox=True, initialization_method="heuristic"
+        ses**2, use_boxcox=True, initialization_method="heuristic"
     )
     with pytest.raises(ValueError, match="use_boxcox was set at model"):
         mod.fit(use_boxcox=True)
@@ -2111,3 +2114,40 @@ def test_invalid_index(reset_randomstate):
     fitted = model.fit(optimized=True, use_brute=True)
     with pytest.warns(ValueWarning, match="No supported"):
         fitted.forecast(steps=157200)
+
+
+def test_initial_level():
+    # GH 8634
+    series = [0.0, 0.0, 0.0, 100.0, 0.0, 0.0, 0.0]
+    es = ExponentialSmoothing(
+        series, initialization_method="known", initial_level=20.0
+    )
+    es_fit = es.fit()
+    es_fit.params
+    assert_allclose(es_fit.params["initial_level"], 20.0)
+
+
+def test_all_initial_values():
+    fit1 = ExponentialSmoothing(
+        aust,
+        seasonal_periods=4,
+        trend="add",
+        seasonal="mul",
+        initialization_method="estimated",
+    ).fit()
+    lvl = np.round(fit1.params["initial_level"])
+    trend = np.round(fit1.params["initial_trend"], 1)
+    seas = np.round(fit1.params["initial_seasons"], 1)
+    fit2 = ExponentialSmoothing(
+        aust,
+        seasonal_periods=4,
+        trend="add",
+        seasonal="mul",
+        initialization_method="known",
+        initial_level=lvl,
+        initial_trend=trend,
+        initial_seasonal=seas,
+    ).fit()
+    assert_allclose(fit2.params["initial_level"], lvl)
+    assert_allclose(fit2.params["initial_trend"], trend)
+    assert_allclose(fit2.params["initial_seasons"], seas)