```python
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
%matplotlib inline
from statsmodels.tsa.exponential_smoothing.ets import ETSModel
```


```python
plt.rcParams['figure.figsize'] = (12, 8)
```

# ETS models

The ETS models are a family of time series models with an underlying state space model consisting of a level component, a trend component (T), a seasonal component (S), and an error term (E).

This notebook gives a very brief introduction to these models and shows how they can be used with statsmodels. For a more thorough treatment we refer to [1], chapter 8 (free online resource), on which the implementation in statsmodels and the examples used in this notebook are based.

[1] Hyndman, Rob J., and George Athanasopoulos. *Forecasting: principles and practice*, 3rd edition, OTexts, 2019. https://www.otexts.org/fpp3/7 

## Simple exponential smoothing

The simplest of the ETS models is also known as *simple exponential smoothing*. In ETS terms, it corresponds to the (A, N, N) model, that is, a model with additive errors, no trend, and no seasonality. The state space formulation of Holt's method is:

\begin{align}
y_{t} &= y_{t-1} + e_t\\
l_{t} &= l_{t-1} + \alpha e_t\\
\end{align}

This state space formulation can be turned into a different formulation, a forecast and a smoothing equation (as can be done with all ETS models):

\begin{align}
\hat{y}_{t|t-1} &= l_{t-1}\\
l_{t} &= \alpha y_{t-1} + (1 - \alpha) l_{t-1}
\end{align}

Here, $\hat{y}_{t|t-1}$ is the forecast/expectation of $y_t$ given the information of the previous step. In the simple exponential smoothing model, the forecast corresponds to the previous level. The second equation (smoothing equation) calculates the next level as weighted average of the previous level and the previous observation.


```python
oildata = [
    111.0091, 130.8284, 141.2871, 154.2278,
    162.7409, 192.1665, 240.7997, 304.2174,
    384.0046, 429.6622, 359.3169, 437.2519,
    468.4008, 424.4353, 487.9794, 509.8284,
    506.3473, 340.1842, 240.2589, 219.0328,
    172.0747, 252.5901, 221.0711, 276.5188,
    271.1480, 342.6186, 428.3558, 442.3946,
    432.7851, 437.2497, 437.2092, 445.3641,
    453.1950, 454.4096, 422.3789, 456.0371,
    440.3866, 425.1944, 486.2052, 500.4291,
    521.2759, 508.9476, 488.8889, 509.8706,
    456.7229, 473.8166, 525.9509, 549.8338,
    542.3405
]
oil = pd.Series(oildata, index=pd.date_range('1965', '2013', freq='AS'))
oil.plot()
plt.ylabel("Annual oil production in Saudi Arabia (Mt)");
```

The plot above shows annual oil production in Saudia Arabia in million tonnes. The data are taken from the R package `fpp2` (companion package to prior version [1]).
Below you can see how to fit a simple exponential smoothing model using statsmodel's ETS implementation to this data. Additionally, the fit using `forecast` in R is shown as comparison.


```python
model = ETSModel(oil, error='add', trend='add', damped_trend=True)
fit = model.fit(maxiter=10000)
oil.plot(label='data')
fit.fittedvalues.plot(label='statsmodels fit')
plt.ylabel("Annual oil production in Saudi Arabia (Mt)");

# obtained from R
params_R = [0.99989969, 0.11888177503085334, 0.80000197, 36.46466837, 34.72584983]
yhat = model.smooth(params_R).fittedvalues
yhat.plot(label='R fit', linestyle='--')

plt.legend();
```

By default the initial states are considered to be fitting parameters and are estimated by maximizing log-likelihood. Additionally it is possible to only use a heuristic for the initial values. In this case this leads to better agreement with the R implementation.


```python
model_heuristic = ETSModel(oil, error='add', trend='add', damped_trend=True,
                          initialization_method='heuristic')
fit_heuristic = model_heuristic.fit()
oil.plot(label='data')
fit.fittedvalues.plot(label='estimated')
fit_heuristic.fittedvalues.plot(label='heuristic', linestyle='--')
plt.ylabel("Annual oil production in Saudi Arabia (Mt)");

# obtained from R
params = [0.99989969, 0.11888177503085334, 0.80000197, 36.46466837, 34.72584983]
yhat = model.smooth(params).fittedvalues
yhat.plot(label='with R params', linestyle=':')

plt.legend();
```


```python
fit.summary()
```

    /home/samuel/Documents/OpenSource/statsmodels/statsmodels/base/model.py:1354: RuntimeWarning: invalid value encountered in sqrt
      bse_ = np.sqrt(np.diag(self.cov_params()))





<table class="simpletable">
<caption>ETS Results</caption>
<tr>
  <th>Dep. Variable:</th>           <td>y</td>        <th>  No. Observations:  </th>    <td>49</td>   
</tr>
<tr>
  <th>Model:</th>               <td>ETS(AAdN)</td>    <th>  Log Likelihood     </th> <td>-258.003</td>
</tr>
<tr>
  <th>Date:</th>            <td>Sun, 07 Jun 2020</td> <th>  AIC                </th>  <td>526.005</td>
</tr>
<tr>
  <th>Time:</th>                <td>19:42:44</td>     <th>  BIC                </th>  <td>535.465</td>
</tr>
<tr>
  <th>Sample:</th>             <td>01-01-1965</td>    <th>  HQIC               </th>  <td>529.594</td>
</tr>
<tr>
  <th></th>                   <td>- 01-01-2013</td>   <th>                     </th>     <td> </td>   
</tr>
<tr>
  <th>Covariance Type:</th>      <td>approx</td>      <th>                     </th>     <td> </td>   
</tr>
</table>
<table class="simpletable">
<tr>
         <td></td>            <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  
</tr>
<tr>
  <th>smoothing_level</th> <td>    1.0000</td> <td>    0.136</td> <td>    7.349</td> <td> 0.000</td> <td>    0.733</td> <td>    1.267</td>
</tr>
<tr>
  <th>smoothing_trend</th> <td>     1e-05</td> <td>      nan</td> <td>      nan</td> <td>   nan</td> <td>      nan</td> <td>      nan</td>
</tr>
<tr>
  <th>damping_trend</th>   <td>    0.9156</td> <td>    0.057</td> <td>   16.115</td> <td> 0.000</td> <td>    0.804</td> <td>    1.027</td>
</tr>
<tr>
  <th>initial_level</th>   <td>   78.0300</td> <td>   52.680</td> <td>    1.481</td> <td> 0.139</td> <td>  -25.220</td> <td>  181.280</td>
</tr>
<tr>
  <th>initial_trend</th>   <td>   36.0141</td> <td>   27.535</td> <td>    1.308</td> <td> 0.191</td> <td>  -17.953</td> <td>   89.981</td>
</tr>
</table>
<table class="simpletable">
<tr>
  <th>Ljung-Box (Q):</th>          <td>nan</td> <th>  Jarque-Bera (JB):  </th> <td>nan</td>
</tr>
<tr>
  <th>Prob(Q):</th>                <td>nan</td> <th>  Prob(JB):          </th> <td>nan</td>
</tr>
<tr>
  <th>Heteroskedasticity (H):</th> <td>nan</td> <th>  Skew:              </th> <td>nan</td>
</tr>
<tr>
  <th>Prob(H) (two-sided):</th>    <td>nan</td> <th>  Kurtosis:          </th> <td>nan</td>
</tr>
</table><br/><br/>Warnings:<br/>[1] Covariance matrix calculated using numerical (complex-step) differentiation.




```python
fit_heuristic.summary()
```




<table class="simpletable">
<caption>ETS Results</caption>
<tr>
  <th>Dep. Variable:</th>           <td>y</td>        <th>  No. Observations:  </th>    <td>49</td>   
</tr>
<tr>
  <th>Model:</th>               <td>ETS(AAdN)</td>    <th>  Log Likelihood     </th> <td>-258.880</td>
</tr>
<tr>
  <th>Date:</th>            <td>Sun, 07 Jun 2020</td> <th>  AIC                </th>  <td>523.760</td>
</tr>
<tr>
  <th>Time:</th>                <td>19:42:44</td>     <th>  BIC                </th>  <td>529.436</td>
</tr>
<tr>
  <th>Sample:</th>             <td>01-01-1965</td>    <th>  HQIC               </th>  <td>525.913</td>
</tr>
<tr>
  <th></th>                   <td>- 01-01-2013</td>   <th>                     </th>     <td> </td>   
</tr>
<tr>
  <th>Covariance Type:</th>      <td>approx</td>      <th>                     </th>     <td> </td>   
</tr>
</table>
<table class="simpletable">
<tr>
         <td></td>            <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  
</tr>
<tr>
  <th>smoothing_level</th> <td>    1.0000</td> <td>    0.097</td> <td>   10.281</td> <td> 0.000</td> <td>    0.809</td> <td>    1.191</td>
</tr>
<tr>
  <th>smoothing_trend</th> <td>    0.1182</td> <td>      nan</td> <td>      nan</td> <td>   nan</td> <td>      nan</td> <td>      nan</td>
</tr>
<tr>
  <th>damping_trend</th>   <td>    0.8000</td> <td>      nan</td> <td>      nan</td> <td>   nan</td> <td>      nan</td> <td>      nan</td>
</tr>
</table>
<table class="simpletable">
<tr>
        <td></td>        <th>initialization method: heuristic</th>
</tr>
<tr>
  <th>initial_level</th>            <td>   33.6309</td>           
</tr>
<tr>
  <th>initial_trend</th>            <td>   34.8115</td>           
</tr>
</table>
<table class="simpletable">
<tr>
  <th>Ljung-Box (Q):</th>          <td>nan</td> <th>  Jarque-Bera (JB):  </th> <td>nan</td>
</tr>
<tr>
  <th>Prob(Q):</th>                <td>nan</td> <th>  Prob(JB):          </th> <td>nan</td>
</tr>
<tr>
  <th>Heteroskedasticity (H):</th> <td>nan</td> <th>  Skew:              </th> <td>nan</td>
</tr>
<tr>
  <th>Prob(H) (two-sided):</th>    <td>nan</td> <th>  Kurtosis:          </th> <td>nan</td>
</tr>
</table><br/><br/>Warnings:<br/>[1] Covariance matrix calculated using numerical (complex-step) differentiation.



## Holt-Winters' seasonal method

The exponential smoothing method can be modified to incorporate a trend and a seasonal component. In the additive Holt-Winters' method, the seasonal component is added to the rest. This model corresponds to the ETS(A, A, A) model, and has the following state space formulation:

\begin{align}
y_t &= l_{t-1} + b_{t-1} + s_{t-m} + e_t\\
l_{t} &= l_{t-1} + b_{t-1} + \alpha e_t\\
b_{t} &= b_{t-1} + \beta e_t\\
s_{t} &= s_{t-m} + \gamma e_t
\end{align}




```python
austourists_data = [
    30.05251300, 19.14849600, 25.31769200, 27.59143700,
    32.07645600, 23.48796100, 28.47594000, 35.12375300, 
    36.83848500, 25.00701700, 30.72223000, 28.69375900, 
    36.64098600, 23.82460900, 29.31168300, 31.77030900,
    35.17787700, 19.77524400, 29.60175000, 34.53884200,
    41.27359900, 26.65586200, 28.27985900, 35.19115300,
    42.20566386, 24.64917133, 32.66733514, 37.25735401,
    45.24246027, 29.35048127, 36.34420728, 41.78208136,
    49.27659843, 31.27540139, 37.85062549, 38.83704413,
    51.23690034, 31.83855162, 41.32342126, 42.79900337,
    55.70835836, 33.40714492, 42.31663797, 45.15712257,
    59.57607996, 34.83733016, 44.84168072, 46.97124960,
    60.01903094, 38.37117851, 46.97586413, 50.73379646,
    61.64687319, 39.29956937, 52.67120908, 54.33231689,
    66.83435838, 40.87118847, 51.82853579, 57.49190993,
    65.25146985, 43.06120822, 54.76075713, 59.83447494,
    73.25702747, 47.69662373, 61.09776802, 66.05576122,
]
index = pd.date_range("1999-03-01", "2015-12-01", freq="3MS")
austourists = pd.Series(austourists_data, index=index)
austourists.plot()
plt.ylabel('Australian Tourists');
```


```python
# fit in statsmodels
model = ETSModel(austourists, error="add", trend="add", seasonal="add",
                damped_trend=True, seasonal_periods=4)
fit = model.fit()

# fit with R params
params_R = [
    0.35445427, 0.03200749, 0.39993387, 0.97999997, 24.01278357, 
    0.97770147, 1.76951063, -0.50735902, -6.61171798, 5.34956637
]
fit_R = model.smooth(params_R)

austourists.plot(label='data')
plt.ylabel('Australian Tourists')

fit.fittedvalues.plot(label='statsmodels fit')
fit_R.fittedvalues.plot(label='R fit', linestyle='--')
plt.legend();
```


```python
fit.summary()
```

    /home/samuel/Documents/OpenSource/statsmodels/statsmodels/base/model.py:1354: RuntimeWarning: invalid value encountered in sqrt
      bse_ = np.sqrt(np.diag(self.cov_params()))





<table class="simpletable">
<caption>ETS Results</caption>
<tr>
  <th>Dep. Variable:</th>           <td>y</td>        <th>  No. Observations:  </th>    <td>68</td>   
</tr>
<tr>
  <th>Model:</th>               <td>ETS(AAdA)</td>    <th>  Log Likelihood     </th> <td>-152.627</td>
</tr>
<tr>
  <th>Date:</th>            <td>Sun, 07 Jun 2020</td> <th>  AIC                </th>  <td>325.254</td>
</tr>
<tr>
  <th>Time:</th>                <td>19:45:10</td>     <th>  BIC                </th>  <td>347.449</td>
</tr>
<tr>
  <th>Sample:</th>             <td>03-01-1999</td>    <th>  HQIC               </th>  <td>334.048</td>
</tr>
<tr>
  <th></th>                   <td>- 12-01-2015</td>   <th>                     </th>     <td> </td>   
</tr>
<tr>
  <th>Covariance Type:</th>      <td>approx</td>      <th>                     </th>     <td> </td>   
</tr>
</table>
<table class="simpletable">
<tr>
           <td></td>             <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  
</tr>
<tr>
  <th>smoothing_level</th>    <td>    0.3398</td> <td>    0.111</td> <td>    3.070</td> <td> 0.002</td> <td>    0.123</td> <td>    0.557</td>
</tr>
<tr>
  <th>smoothing_trend</th>    <td>    0.0259</td> <td>    0.008</td> <td>    3.157</td> <td> 0.002</td> <td>    0.010</td> <td>    0.042</td>
</tr>
<tr>
  <th>smoothing_seasonal</th> <td>    0.4011</td> <td>    0.080</td> <td>    5.041</td> <td> 0.000</td> <td>    0.245</td> <td>    0.557</td>
</tr>
<tr>
  <th>damping_trend</th>      <td>    0.9800</td> <td>      nan</td> <td>      nan</td> <td>   nan</td> <td>      nan</td> <td>      nan</td>
</tr>
<tr>
  <th>initial_level</th>      <td>   24.6362</td> <td>      nan</td> <td>      nan</td> <td>   nan</td> <td>      nan</td> <td>      nan</td>
</tr>
<tr>
  <th>initial_trend</th>      <td>    0.6148</td> <td>    0.392</td> <td>    1.568</td> <td> 0.117</td> <td>   -0.154</td> <td>    1.383</td>
</tr>
<tr>
  <th>initial_seasonal.0</th> <td>    1.3700</td> <td>      nan</td> <td>      nan</td> <td>   nan</td> <td>      nan</td> <td>      nan</td>
</tr>
<tr>
  <th>initial_seasonal.1</th> <td>   -1.1431</td> <td>      nan</td> <td>      nan</td> <td>   nan</td> <td>      nan</td> <td>      nan</td>
</tr>
<tr>
  <th>initial_seasonal.2</th> <td>   -6.6725</td> <td>      nan</td> <td>      nan</td> <td>   nan</td> <td>      nan</td> <td>      nan</td>
</tr>
<tr>
  <th>initial_seasonal.3</th> <td>    4.8117</td> <td>      nan</td> <td>      nan</td> <td>   nan</td> <td>      nan</td> <td>      nan</td>
</tr>
</table>
<table class="simpletable">
<tr>
  <th>Ljung-Box (Q):</th>          <td>nan</td> <th>  Jarque-Bera (JB):  </th> <td>nan</td>
</tr>
<tr>
  <th>Prob(Q):</th>                <td>nan</td> <th>  Prob(JB):          </th> <td>nan</td>
</tr>
<tr>
  <th>Heteroskedasticity (H):</th> <td>nan</td> <th>  Skew:              </th> <td>nan</td>
</tr>
<tr>
  <th>Prob(H) (two-sided):</th>    <td>nan</td> <th>  Kurtosis:          </th> <td>nan</td>
</tr>
</table><br/><br/>Warnings:<br/>[1] Covariance matrix calculated using numerical (complex-step) differentiation.




```python
fit._rank
```




    10




```python

```

e4091700f7903ecd6a8e1e4089ec0baeb38f3d0d
ipynb
Jupyter Notebook
examples/notebooks/ets.ipynb
ronenmr/statsmodels
7de865a99388047d6dc37ef8a06e943bd8211a00
BSD-3-Clause
examples/notebooks/ets.ipynb
ronenmr/statsmodels
7de865a99388047d6dc37ef8a06e943bd8211a00
BSD-3-Clause
examples/notebooks/ets.ipynb
ronenmr/statsmodels
7de865a99388047d6dc37ef8a06e943bd8211a00
BSD-3-Clause
Qwen/Qwen-72B
1. YES
2. YES