# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: laughingrook
#     language: python
#     name: laughingrook
# ---

# # Holt's Exponential Smoothing
# > Use Holt's exponential smoothing method to model a seasonal time series with no trend.
#
# - toc: true
# - badges: true
# - comments: true
# - categories: [m249, timeseries]
# - permalink: /2022/06/06/holts_exponential_smoothing_time_series/

# ## Dependencies

import numpy as np
import pandas as pd
from scipy import stats as st
from statsmodels import api as sm
from statsmodels import datasets
from statsmodels.tsa import api as tsa
from statsmodels.graphics import tsaplots
from matplotlib import pyplot as plt
import seaborn as sns

# Set the graphing defaults.
# (This is optional.)

sns.set_theme()
plt.rc("figure", figsize=(8, 6))
plt.rc("font", size=14)

# ## Notes
#
# - Data was sourced from the Rdatasets
#    - <NAME> (2022). Rdatasets: A collection of datasets originally distributed in various R packages. R package version 1.0.0, https://vincentarelbundock.github.io/Rdatasets.
# - `Dataset` is a class from StatsModels
#    - See [The Datasets Package](https://www.statsmodels.org/devel/datasets/index.html)
# - Time series is seasonal and has no marked updward or downward trend, so it is reasonable to use Holt's exponential smoothing
# - Model parameters
#    - $\alpha \simeq$ 0.985050 $\to$ more weight is placed on more recent values of the *level*
#    - $\gamma \simeq$ 0.418600 $\to$ approximately equal weighting placed on distant past and recent values of the *slope*
# - Residuals
#    - No sample autocorrelations exceed the 95% signficance boundaries
#    - Ljung-Box test yielded a test statistic = 24.420038, *p* = 0.225
#       - Little to no evidence against a white noise model
#    - Time plot of forecasting errors shows no systematic change in level or variance
#    - Histogram of residuals shows assumption that the forecast errors are normally distributed appears plausible

# ## Inputs/Outputs
#
# ### Inputs
#
# - *hartnagel* (`Dataset`), Canadian Crime-Rates Time Series
#
# ### Intermediate
#
# - *ts* (`Series`), annual total fertility rate per 1000 women with `DatetimeIndex`
# - *holt* (`Holt`), input *ts* modelled using Holt's exponential smoothing
#    - See [Holt](https://www.statsmodels.org/dev/generated/statsmodels.tsa.holtwinters.Holt.html)
#
# ### Outputs
#
# - *fitted* (`HoltWintersResults`), modelled time series
#    - See [HoltWintersResults](https://www.statsmodels.org/dev/generated/statsmodels.tsa.holtwinters.HoltWintersResults.html)
# - *fcast* (`Series`), quantified forecast for 1969

# ## Initialise the time series

# ### Load the data

hartnagel = datasets.get_rdataset('Hartnagel', package='carData', cache=True)
hartnagel.data.info()

# ### Initialise the Series (Pandas)
#
# The first observation is in 1931.

hartnagel.data['year'].head(1)

ts = pd.Series(
    data=hartnagel.data['tfr'].to_numpy(),
    name='TOTAL FERTILITY RATE',
    index=pd.date_range(
        start='1931',
        periods=hartnagel.data['tfr'].size,
        freq='A-DEC',
        name='YEAR'
    )
)
ts.info()

# ### Plot the time series

sns.lineplot(x=ts.index, y=ts, lw=2, label='obs')
plt.ylabel(ts.name)
plt.legend()
plt.show()

# ## Smooth the time series

holt = tsa.Holt(ts, initialization_method='estimated')
fitted = holt.fit()

# Plot the observerved and fitted time series.

sns.lineplot(x=ts.index, y=ts, lw=2, label='obs')
fitted.fittedvalues.plot(lw=2, label='fit')
plt.ylabel(ts.name)
plt.legend()
plt.show()

# ### Model parameters

params = pd.Series(
    data={
        'alpha': fitted.params['smoothing_level'],
        'gamma': fitted.params['smoothing_trend'],
        'sse': fitted.sse,
    },
    name='smoothing parameters'
)
params

# ### Forecasting errors

# #### Time series of the forecasting errors

sns.lineplot(x=ts.index, y=fitted.resid, lw=2)
plt.ylabel('RESIDUAL')
plt.show()

# #### In-sample autocorrelations

tsaplots.plot_acf(fitted.resid, lags=20, zero=False)
plt.xticks(range(21))
plt.show()

# #### Ljung-Box test

sm.stats.acorr_ljungbox(fitted.resid, lags=[20])

# #### Histogram of the forecasting errors

sns.histplot(x=fitted.resid, kde=True)
plt.show()

# ### Forecasting

# #### Plot of forecasts (+3 years)

fcast3 = fitted.forecast(3)
sns.lineplot(x=ts.index, y=ts, lw=2, label='obs')
fitted.fittedvalues.plot(lw=2, label='fit')
fcast3.plot(lw=2, label='fcast')
plt.axvline(fcast3.index[0], alpha=0.7, lw=0.5, ls='--', color='black')
plt.legend()
plt.ylabel(ts.name)
plt.show()

# #### Quantify a forecast

fcast = fitted.forecast(1)[0]
std_err = np.sqrt(
    st.norm().ppf(0.975) * (fitted.sse / fitted.fittedvalues.size)
)
pd.Series(
    data={'pred': fcast, 'lcb': fcast - std_err, 'ucb': fcast + std_err},
    name='1969 predicted fertility rate'
)

# ## References
#
# - M249, Book 2, Part 2.7 & 2.9

StarcoderdataJupyterScriptsDedupFiltered