Repository: statsmodels/statsmodels
Pull request: ENH: added diagnostics test to ETS model
- [ ] closes #6889
- [ ] tests added / passed. 
- [ ] code/documentation is well formatted.  
- [ ] properly formatted commit message. See 
      [NumPy's guide](https://docs.scipy.org/doc/numpy-1.15.1/dev/gitwash/development_workflow.html#writing-the-commit-message). 

I added the diagnostic test methods to the base class. In a comparison with the statespace implementation the heteroskedasticity is still different, I am unsure why (test fails currently, if this behaviour is expected I can adapt the test).

Additionally, I modified initial seasonal values as described in #6889.

Problem statement:
ENH: Remove NaNs from ETS summary
When fitting the AAdA model to the austourists data there's still some NaNs showing up in the summary frame.

- The final section is all NaNs, because I didn't implement these yet for ETS.
- The error of `damping_trend` is NaN, because it is directly on the fitting bounds
- The error of `initial_level` and the `initial_seasonal`s is NaN. Is this a bug? I saw that the statespace implementations uses only 3 initial seasonal values - so I should probably also fix one of the initial seasonal values to a constant value? I thought about fixing the last initial seasonal value to 0 ( or 1), subtracting it from (dividing by) the other initial seasonal values, and adding this value to the initial level (multiplying the initial level by this value) for the additive (multiplicative) models. @ChadFulton, do you have some ideas/recommendations?

Summary Frame:
<details>
```
                                 ETS Results                                  
==============================================================================
Dep. Variable:                      y   No. Observations:                   68
Model:                      ETS(AAdA)   Log Likelihood                -152.627
Date:                Sat, 18 Jul 2020   AIC                            325.254
Time:                        10:29:54   BIC                            347.449
Sample:                    03-01-1999   HQIC                           334.048
                         - 12-01-2015   Scale                            5.213
Covariance Type:               approx                                         
======================================================================================
                         coef    std err          z      P>|z|      [0.025      0.975]
--------------------------------------------------------------------------------------
smoothing_level        0.3399      0.111      3.070      0.002       0.123       0.557
smoothing_trend        0.0259      0.008      3.156      0.002       0.010       0.042
smoothing_seasonal     0.4011      0.080      5.041      0.000       0.245       0.557
damping_trend          0.9800        nan        nan        nan         nan         nan
initial_level         24.6361        nan        nan        nan         nan         nan
initial_trend          0.6147      0.392      1.567      0.117      -0.154       1.383
initial_seasonal.0     1.3698        nan        nan        nan         nan         nan
initial_seasonal.1    -1.1431        nan        nan        nan         nan         nan
initial_seasonal.2    -6.6727        nan        nan        nan         nan         nan
initial_seasonal.3     4.8119        nan        nan        nan         nan         nan
===================================================================================
Ljung-Box (Q):                         nan   Jarque-Bera (JB):                  nan
Prob(Q):                               nan   Prob(JB):                          nan
Heteroskedasticity (H):                nan   Skew:                              nan
Prob(H) (two-sided):                   nan   Kurtosis:                          nan
===================================================================================

Warnings:
[1] Covariance matrix calculated using numerical (complex-step) differentiation.
```
</details>

Hints from issue discussion:

Commit message: fix merge

Patch to review:
diff --git a/statsmodels/tsa/exponential_smoothing/base.py b/statsmodels/tsa/exponential_smoothing/base.py
index 7c7e64b076e..3ac4c5c6bb1 100644
--- a/statsmodels/tsa/exponential_smoothing/base.py
+++ b/statsmodels/tsa/exponential_smoothing/base.py
@@ -347,10 +347,6 @@ def __init__(self, model, params, scale=1.0):
     def nobs_effective(self):
         raise NotImplementedError
 
-    @cache_readonly
-    def df_model(self):
-        raise NotImplementedError
-
     @cache_readonly
     def df_resid(self):
         return self.nobs_effective - self.df_model
@@ -501,6 +497,333 @@ def cov_params_approx(self):
             self._cov_approx_complex_step, self._cov_approx_centered
         )
 
+    def test_serial_correlation(self, method, lags=None):
+        """
+        Ljung-Box test for no serial correlation of standardized residuals
+
+        Null hypothesis is no serial correlation.
+
+        Parameters
+        ----------
+        method : {'ljungbox','boxpierece', None}
+            The statistical test for serial correlation. If None, an attempt is
+            made to select an appropriate test.
+        lags : None, int or array_like
+            If lags is an integer then this is taken to be the largest lag
+            that is included, the test result is reported for all smaller lag
+            length.
+            If lags is a list or array, then all lags are included up to the
+            largest lag in the list, however only the tests for the lags in the
+            list are reported.
+            If lags is None, then the default maxlag is 12*(nobs/100)^{1/4}
+
+        Returns
+        -------
+        output : ndarray
+            An array with `(test_statistic, pvalue)` for each endogenous
+            variable and each lag. The array is then sized
+            `(k_endog, 2, lags)`. If the method is called as
+            `ljungbox = res.test_serial_correlation()`, then `ljungbox[i]`
+            holds the results of the Ljung-Box test (as would be returned by
+            `statsmodels.stats.diagnostic.acorr_ljungbox`) for the `i` th
+            endogenous variable.
+
+        See Also
+        --------
+        statsmodels.stats.diagnostic.acorr_ljungbox
+            Ljung-Box test for serial correlation.
+
+        Notes
+        -----
+        For statespace models: let `d` = max(loglikelihood_burn, nobs_diffuse);
+        this test is calculated ignoring the first `d` residuals.
+
+        Output is nan for any endogenous variable which has missing values.
+        """
+        if method is None:
+            method = 'ljungbox'
+
+        if self.standardized_forecasts_error is None:
+            raise ValueError('Cannot compute test statistic when standardized'
+                             ' forecast errors have not been computed.')
+
+        if method == 'ljungbox' or method == 'boxpierce':
+            from statsmodels.stats.diagnostic import acorr_ljungbox
+            if hasattr(self, "loglikelihood_burn"):
+                d = np.maximum(self.loglikelihood_burn, self.nobs_diffuse)
+                # This differs from self.nobs_effective because here we want to
+                # exclude exact diffuse periods, whereas self.nobs_effective
+                # only excludes explicitly burned (usually approximate diffuse)
+                # periods.
+                nobs_effective = self.nobs - d
+            else:
+                nobs_effective = self.nobs_effective
+            output = []
+
+            # Default lags for acorr_ljungbox is 40, but may not always have
+            # that many observations
+            if lags is None:
+                seasonal_periods = getattr(self.model, "seasonal_periods", 0)
+                if seasonal_periods:
+                    lags = min(2 * seasonal_periods, nobs_effective // 5)
+                else:
+                    lags = min(10, nobs_effective // 5)
+
+                warnings.warn(
+                    "The default value of lags is changing.  After 0.12, "
+                    "this value will become min(10, nobs//5) for non-seasonal "
+                    "time series and min (2*m, nobs//5) for seasonal time "
+                    "series. Directly set lags to silence this warning.",
+                    FutureWarning
+                )
+
+            for i in range(self.model.k_endog):
+                if hasattr(self, "filter_results"):
+                    x = self.filter_results.standardized_forecasts_error[i][d:]
+                else:
+                    x = self.standardized_forecasts_error
+                results = acorr_ljungbox(
+                    x, lags=lags, boxpierce=(method == 'boxpierce'),
+                    return_df=False)
+                if method == 'ljungbox':
+                    output.append(results[0:2])
+                else:
+                    output.append(results[2:])
+
+            output = np.c_[output]
+        else:
+            raise NotImplementedError('Invalid serial correlation test'
+                                      ' method.')
+        return output
+
+    def test_heteroskedasticity(self, method, alternative='two-sided',
+                                use_f=True):
+        r"""
+        Test for heteroskedasticity of standardized residuals
+
+        Tests whether the sum-of-squares in the first third of the sample is
+        significantly different than the sum-of-squares in the last third
+        of the sample. Analogous to a Goldfeld-Quandt test. The null hypothesis
+        is of no heteroskedasticity.
+
+        Parameters
+        ----------
+        method : {'breakvar', None}
+            The statistical test for heteroskedasticity. Must be 'breakvar'
+            for test of a break in the variance. If None, an attempt is
+            made to select an appropriate test.
+        alternative : str, 'increasing', 'decreasing' or 'two-sided'
+            This specifies the alternative for the p-value calculation. Default
+            is two-sided.
+        use_f : bool, optional
+            Whether or not to compare against the asymptotic distribution
+            (chi-squared) or the approximate small-sample distribution (F).
+            Default is True (i.e. default is to compare against an F
+            distribution).
+
+        Returns
+        -------
+        output : ndarray
+            An array with `(test_statistic, pvalue)` for each endogenous
+            variable. The array is then sized `(k_endog, 2)`. If the method is
+            called as `het = res.test_heteroskedasticity()`, then `het[0]` is
+            an array of size 2 corresponding to the first endogenous variable,
+            where `het[0][0]` is the test statistic, and `het[0][1]` is the
+            p-value.
+
+        Notes
+        -----
+        The null hypothesis is of no heteroskedasticity. That means different
+        things depending on which alternative is selected:
+
+        - Increasing: Null hypothesis is that the variance is not increasing
+          throughout the sample; that the sum-of-squares in the later
+          subsample is *not* greater than the sum-of-squares in the earlier
+          subsample.
+        - Decreasing: Null hypothesis is that the variance is not decreasing
+          throughout the sample; that the sum-of-squares in the earlier
+          subsample is *not* greater than the sum-of-squares in the later
+          subsample.
+        - Two-sided: Null hypothesis is that the variance is not changing
+          throughout the sample. Both that the sum-of-squares in the earlier
+          subsample is not greater than the sum-of-squares in the later
+          subsample *and* that the sum-of-squares in the later subsample is
+          not greater than the sum-of-squares in the earlier subsample.
+
+        For :math:`h = [T/3]`, the test statistic is:
+
+        .. math::
+
+            H(h) = \sum_{t=T-h+1}^T  \tilde v_t^2
+            \Bigg / \sum_{t=d+1}^{d+1+h} \tilde v_t^2
+
+        where :math:`d` = max(loglikelihood_burn, nobs_diffuse)` (usually
+        corresponding to diffuse initialization under either the approximate
+        or exact approach).
+
+        This statistic can be tested against an :math:`F(h,h)` distribution.
+        Alternatively, :math:`h H(h)` is asymptotically distributed according
+        to :math:`\chi_h^2`; this second test can be applied by passing
+        `asymptotic=True` as an argument.
+
+        See section 5.4 of [1]_ for the above formula and discussion, as well
+        as additional details.
+
+        TODO
+
+        - Allow specification of :math:`h`
+
+        References
+        ----------
+        .. [1] Harvey, Andrew C. 1990. *Forecasting, Structural Time Series*
+               *Models and the Kalman Filter.* Cambridge University Press.
+        """
+        if method is None:
+            method = 'breakvar'
+
+        if self.standardized_forecasts_error is None:
+            raise ValueError('Cannot compute test statistic when standardized'
+                             ' forecast errors have not been computed.')
+
+        if method == 'breakvar':
+            # Store some values
+            if hasattr(self, "filter_results"):
+                squared_resid = (
+                    self.filter_results.standardized_forecasts_error**2
+                )
+                d = np.maximum(self.loglikelihood_burn, self.nobs_diffuse)
+                # This differs from self.nobs_effective because here we want to
+                # exclude exact diffuse periods, whereas self.nobs_effective
+                # only excludes explicitly burned (usually approximate diffuse)
+                # periods.
+                nobs_effective = self.nobs - d
+            else:
+                squared_resid = self.standardized_forecasts_error**2
+                if squared_resid.ndim == 1:
+                    squared_resid = squared_resid[np.newaxis, :]
+                nobs_effective = self.nobs_effective
+                d = 0
+
+            test_statistics = []
+            p_values = []
+            for i in range(self.model.k_endog):
+                h = int(np.round(nobs_effective / 3))
+                numer_resid = squared_resid[i, -h:]
+                numer_resid = numer_resid[~np.isnan(numer_resid)]
+                numer_dof = len(numer_resid)
+
+                denom_resid = squared_resid[i, d:d + h]
+                denom_resid = denom_resid[~np.isnan(denom_resid)]
+                denom_dof = len(denom_resid)
+
+                if numer_dof < 2:
+                    warnings.warn('Early subset of data for variable %d'
+                                  '  has too few non-missing observations to'
+                                  ' calculate test statistic.' % i)
+                    numer_resid = np.nan
+                if denom_dof < 2:
+                    warnings.warn('Later subset of data for variable %d'
+                                  '  has too few non-missing observations to'
+                                  ' calculate test statistic.' % i)
+                    denom_resid = np.nan
+
+                test_statistic = np.sum(numer_resid) / np.sum(denom_resid)
+
+                # Setup functions to calculate the p-values
+                if use_f:
+                    from scipy.stats import f
+                    pval_lower = lambda test_statistics: f.cdf(  # noqa:E731
+                        test_statistics, numer_dof, denom_dof)
+                    pval_upper = lambda test_statistics: f.sf(  # noqa:E731
+                        test_statistics, numer_dof, denom_dof)
+                else:
+                    from scipy.stats import chi2
+                    pval_lower = lambda test_statistics: chi2.cdf(  # noqa:E731
+                        numer_dof * test_statistics, denom_dof)
+                    pval_upper = lambda test_statistics: chi2.sf(  # noqa:E731
+                        numer_dof * test_statistics, denom_dof)
+
+                # Calculate the one- or two-sided p-values
+                alternative = alternative.lower()
+                if alternative in ['i', 'inc', 'increasing']:
+                    p_value = pval_upper(test_statistic)
+                elif alternative in ['d', 'dec', 'decreasing']:
+                    test_statistic = 1. / test_statistic
+                    p_value = pval_upper(test_statistic)
+                elif alternative in ['2', '2-sided', 'two-sided']:
+                    p_value = 2 * np.minimum(
+                        pval_lower(test_statistic),
+                        pval_upper(test_statistic)
+                    )
+                else:
+                    raise ValueError('Invalid alternative.')
+
+                test_statistics.append(test_statistic)
+                p_values.append(p_value)
+
+            output = np.c_[test_statistics, p_values]
+        else:
+            raise NotImplementedError('Invalid heteroskedasticity test'
+                                      ' method.')
+
+        return output
+
+    def test_normality(self, method):
+        """
+        Test for normality of standardized residuals.
+
+        Null hypothesis is normality.
+
+        Parameters
+        ----------
+        method : {'jarquebera', None}
+            The statistical test for normality. Must be 'jarquebera' for
+            Jarque-Bera normality test. If None, an attempt is made to select
+            an appropriate test.
+
+        See Also
+        --------
+        statsmodels.stats.stattools.jarque_bera
+            The Jarque-Bera test of normality.
+
+        Notes
+        -----
+        For statespace models: let `d` = max(loglikelihood_burn, nobs_diffuse);
+        this test is calculated ignoring the first `d` residuals.
+
+        In the case of missing data, the maintained hypothesis is that the
+        data are missing completely at random. This test is then run on the
+        standardized residuals excluding those corresponding to missing
+        observations.
+        """
+        if method is None:
+            method = 'jarquebera'
+
+        if self.standardized_forecasts_error is None:
+            raise ValueError('Cannot compute test statistic when standardized'
+                             ' forecast errors have not been computed.')
+
+        if method == 'jarquebera':
+            from statsmodels.stats.stattools import jarque_bera
+            if hasattr(self, "loglikelihood_burn"):
+                d = np.maximum(self.loglikelihood_burn, self.nobs_diffuse)
+            else:
+                d = 0
+            output = []
+            for i in range(self.model.k_endog):
+                if hasattr(self, "fiter_results"):
+                    resid = self.filter_results.standardized_forecasts_error[
+                        i, d:
+                    ]
+                else:
+                    resid = self.standardized_forecasts_error
+                mask = ~np.isnan(resid)
+                output.append(jarque_bera(resid[mask]))
+        else:
+            raise NotImplementedError('Invalid normality test method.')
+
+        return np.array(output)
+
     def summary(
         self,
         alpha=0.05,
diff --git a/statsmodels/tsa/exponential_smoothing/ets.py b/statsmodels/tsa/exponential_smoothing/ets.py
index aec58ba6dfa..8abd51abfd7 100644
--- a/statsmodels/tsa/exponential_smoothing/ets.py
+++ b/statsmodels/tsa/exponential_smoothing/ets.py
@@ -195,6 +195,14 @@
 #   - set up the constraint bounds and function
 # * Since the traditional bounds are nonlinear for beta and gamma, if no bounds
 #   are given, we internally use beta_star and gamma_star for fitting
+# * When estimating initial level and initial seasonal values, one of them has
+#   to be removed in order to have a well posed problem. I am solving this by
+#   fixing the last initial seasonal value to 0 (for additive seasonality) or 1
+#   (for multiplicative seasonality).
+#   For the additive models, this means I have to subtract the last initial
+#   seasonal value from all initial seasonal values and add it to the initial
+#   level; for the multiplicative models I do the same with division and
+#   multiplication
 
 
 class ETSModel(base.StateSpaceMLEModel):
@@ -238,8 +246,8 @@ class ETSModel(base.StateSpaceMLEModel):
     initial_trend : float, optional
         The initial trend component. Only used if initialization is 'known'.
     initial_seasonal : array_like, optional
-        The initial seasonal component. An array of length `seasonal`.  Only
-        used if initialization is 'known'.
+        The initial seasonal component. An array of length `seasonal_periods`.
+        Only used if initialization is 'known'.
     bounds : dict or None, optional
         A dictionary with parameter names as keys and the respective bounds
         intervals as values (lists/tuples/arrays).
@@ -527,8 +535,8 @@ def set_initialization_method(
             The initial trend component. Only used if initialization is
             'known'.
         initial_seasonal : array_like, optional
-            The initial seasonal component. An array of length `seasonal`. Only
-            used if initialization is 'known'.
+            The initial seasonal component. An array of length
+            `seasonal_periods`. Only used if initialization is 'known'.
         """
         self.initialization_method = string_like(
             initialization_method,
@@ -631,6 +639,10 @@ def prepare_data(data):
     def nobs_effective(self):
         return self.nobs
 
+    @property
+    def k_endog(self):
+        return 1
+
     @property
     def short_name(self):
         name = "".join(
@@ -813,7 +825,7 @@ def _default_start_params(self):
     @property
     def _start_params(self):
         """
-        Returns default start params in the format of external parameters.
+        Default start params in the format of external parameters.
         This should not be called directly, but by calling
         ``self.start_params``.
         """
@@ -823,11 +835,21 @@ def _start_params(self):
                 params.append(self._default_start_params[p])
 
         if self.initialization_method == "estimated":
+            lvl_idx = len(params)
             params += [self.initial_level]
             if self.has_trend:
                 params += [self.initial_trend]
             if self.has_seasonal:
-                params += self.initial_seasonal.tolist()
+                # we have to adapt the seasonal values a bit to make sure the
+                # problem is well posed (see implementation notes above)
+                initial_seasonal = self.initial_seasonal
+                if self.seasonal == "mul":
+                    params[lvl_idx] *= initial_seasonal[-1]
+                    initial_seasonal /= initial_seasonal[-1]
+                else:
+                    params[lvl_idx] += initial_seasonal[-1]
+                    initial_seasonal -= initial_seasonal[-1]
+                params += initial_seasonal.tolist()
 
         return np.array(params)
 
@@ -877,8 +899,14 @@ def _setup_bounds(self):
         # other bounds for phi and initial states
         lb[3], ub[3] = 0.8, 0.98
         if self.initialization_method == "estimated":
-            lb[4:] = -np.inf
-            ub[4:] = np.inf
+            lb[4:-1] = -np.inf
+            ub[4:-1] = np.inf
+            # fix the last initial_seasonal to 0 or 1, otherwise the equation
+            # is underdetermined
+            if self.seasonal == "mul":
+                lb[-1], ub[-1] = 1, 1
+            else:
+                lb[-1], ub[-1] = 0, 0
 
         # set lb and ub for parameters with bounds
         for p in self._internal_param_names:
@@ -1305,6 +1333,8 @@ def __init__(self, model, params, results):
         self._llf = model.loglike(params)
         self._residuals = model._residuals(yhat)
         self._fittedvalues = yhat
+        # scale is concentrated in this model formulation and corresponds to
+        # mean squared residuals, see docstring of model.loglike
         scale = np.mean(self._residuals ** 2)
         super().__init__(model, params, scale=scale)
 
@@ -1359,6 +1389,17 @@ def __init__(self, model, params, results):
             self.phi = internal_params[3]
             self.damping_trend = self.phi
 
+        # degrees of freedom of model
+        k_free_params = self.k_params - len(self.fixed_params)
+        self.df_model = k_free_params + 1
+
+        # standardized forecasting error
+        self.mean_resid = np.mean(self.resid)
+        self.scale_resid = np.std(self.resid, ddof=1)
+        self.standardized_forecasts_error = (
+            (self.resid - self.mean_resid) / self.scale_resid
+        )
+
         # Setup covariance matrix notes dictionary
         # For now, only support "approx"
         if not hasattr(self, "cov_kwds"):
@@ -1396,10 +1437,6 @@ def __init__(self, model, params, results):
     def nobs_effective(self):
         return self.nobs
 
-    @cache_readonly
-    def df_model(self):
-        return self.model.k_params
-
     @cache_readonly
     def fittedvalues(self):
         return self._fittedvalues
@@ -1721,8 +1758,7 @@ def simulate(
         y = np.empty((nsimulations, repetitions))
 
         # get random error eps
-        sigma = np.sqrt(self.sse / self.df_resid)
-        # sigma = np.sqrt(self.scale)
+        sigma = np.sqrt(self.scale)
         if isinstance(random_errors, np.ndarray):
             if random_errors.shape != (nsimulations, repetitions):
                 raise ValueError(
diff --git a/statsmodels/tsa/tests/test_exponential_smoothing.py b/statsmodels/tsa/tests/test_exponential_smoothing.py
index 5440b297903..161ea3fbcad 100644
--- a/statsmodels/tsa/tests/test_exponential_smoothing.py
+++ b/statsmodels/tsa/tests/test_exponential_smoothing.py
@@ -719,6 +719,20 @@ def test_results_vs_statespace(statespace_comparison):
         statespace_results.fittedvalues.values
     )
 
+    # compare diagnostics
+    assert_almost_equal(
+        ets_results.test_serial_correlation(method="ljungbox"),
+        statespace_results.test_serial_correlation(method="ljungbox"),
+    )
+    assert_almost_equal(
+        ets_results.test_normality(method="jarquebera"),
+        statespace_results.test_normality(method="jarquebera"),
+    )
+    assert_almost_equal(
+        ets_results.test_heteroskedasticity(method="breakvar"),
+        statespace_results.test_heteroskedasticity(method="breakvar"),
+    )
+
 
 def test_prediction_results_vs_statespace(statespace_comparison):
     ets_results, statespace_results = statespace_comparison


Code review comments:
diff_hunk: @@ -719,6 +719,20 @@ def test_results_vs_statespace(statespace_comparison):
         statespace_results.fittedvalues.values
     )
 
+    # compare diagnostics
+    assert_almost_equal(
path: statsmodels/tsa/tests/test_exponential_smoothing.py
text: @user1:
You need to wrap the assert almost_equal with

```

with pytest.warns(FutureWarning):
   assert_almost_equal(...
```

since this warning is transformed into an error to make sure we aren't leaking any into a test run.