Introduction and Preprocessing

To run these notebooks you will need the sktime library. You can install it with pip or conda

pip install sktime

or

conda install conda-forge::sktime

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sktime.datasets import load_airline, load_uschange, load_basic_motions
from sktime.utils.plotting import plot_series

Intro

Time Series Data

univariate_ts = pd.read_csv("data/airline.csv", index_col=0)
univariate_ts.plot()

<Axes: xlabel='Period'>

_, multivariate_ts = load_uschange()
sns.lineplot(data=multivariate_ts.iloc[:, :2], color="C0")

<Axes: >

multivariate_ts_panel, labels = load_basic_motions()
colors = ["C0", "C1"]
linestyles = ["-", "--", ":", "-."]
for i in range(2):
    for j in range(2):
        multivariate_ts_panel.iloc[j, i].plot(color=colors[i], linestyle=linestyles[i + j], label=f"channel {j+1} of ts {i+1}")
plt.legend()
plt.show()

Load the data

# from csv
ts = pd.read_csv("data/airline.csv", index_col=0)
ts.head()

	Number of airline passengers
Period
1949-01	112.0
1949-02	118.0
1949-03	132.0
1949-04	129.0
1949-05	121.0

# directly from sktime
ts = load_airline()
ts.head()

Period
1949-01    112.0
1949-02    118.0
1949-03    132.0
1949-04    129.0
1949-05    121.0
Freq: M, Name: Number of airline passengers, dtype: float64

Some basic statistics

ts.shape

(144,)

len(ts)

144

ts.index

PeriodIndex(['1949-01', '1949-02', '1949-03', '1949-04', '1949-05', '1949-06',
             '1949-07', '1949-08', '1949-09', '1949-10',
             ...
             '1960-03', '1960-04', '1960-05', '1960-06', '1960-07', '1960-08',
             '1960-09', '1960-10', '1960-11', '1960-12'],
            dtype='period[M]', name='Period', length=144)

# what is the frequency of the time series?
ts.index.freqstr

'M'

ts.index.year

Index([1949, 1949, 1949, 1949, 1949, 1949, 1949, 1949, 1949, 1949,
       ...
       1960, 1960, 1960, 1960, 1960, 1960, 1960, 1960, 1960, 1960],
      dtype='int64', name='Period', length=144)

# how many values from each year?
ts.index.year.value_counts()

Period
1949    12
1950    12
1951    12
1952    12
1953    12
1954    12
1955    12
1956    12
1957    12
1958    12
1959    12
1960    12
Name: count, dtype: int64

ts.array

<NumpyExtensionArray>
[112.0, 118.0, 132.0, 129.0, 121.0, 135.0, 148.0, 148.0, 136.0, 119.0,
 ...
 419.0, 461.0, 472.0, 535.0, 622.0, 606.0, 508.0, 461.0, 390.0, 432.0]
Length: 144, dtype: float64

ts.to_numpy()

array([112., 118., 132., 129., 121., 135., 148., 148., 136., 119., 104.,
       118., 115., 126., 141., 135., 125., 149., 170., 170., 158., 133.,
       114., 140., 145., 150., 178., 163., 172., 178., 199., 199., 184.,
       162., 146., 166., 171., 180., 193., 181., 183., 218., 230., 242.,
       209., 191., 172., 194., 196., 196., 236., 235., 229., 243., 264.,
       272., 237., 211., 180., 201., 204., 188., 235., 227., 234., 264.,
       302., 293., 259., 229., 203., 229., 242., 233., 267., 269., 270.,
       315., 364., 347., 312., 274., 237., 278., 284., 277., 317., 313.,
       318., 374., 413., 405., 355., 306., 271., 306., 315., 301., 356.,
       348., 355., 422., 465., 467., 404., 347., 305., 336., 340., 318.,
       362., 348., 363., 435., 491., 505., 404., 359., 310., 337., 360.,
       342., 406., 396., 420., 472., 548., 559., 463., 407., 362., 405.,
       417., 391., 419., 461., 472., 535., 622., 606., 508., 461., 390.,
       432.])

ts.describe()

count    144.000000
mean     280.298611
std      119.966317
min      104.000000
25%      180.000000
50%      265.500000
75%      360.500000
max      622.000000
Name: Number of airline passengers, dtype: float64

ts.info()

<class 'pandas.core.series.Series'>
PeriodIndex: 144 entries, 1949-01 to 1960-12
Freq: M
Series name: Number of airline passengers
Non-Null Count  Dtype  
--------------  -----  
144 non-null    float64
dtypes: float64(1)
memory usage: 2.2 KB

Plotting

# pandas
ts.plot()
plt.show()

# matplotlib
plt.plot(ts.index.astype("str"), ts.to_numpy())
plt.xticks(ts.index.astype("str")[::20], rotation=45)
plt.show()

# seaborn
sns.lineplot(x=ts.index.astype("str"), y=ts)
plt.xticks(ts.index.astype("str")[::20], rotation=45)
plt.show()

# sktime
plot_series(ts)
plt.show()

/Users/francesco/miniforge3/envs/timeseries_dl/lib/python3.12/site-packages/gluonts/json.py:102: UserWarning: Using `json`-module for json-handling. Consider installing one of `orjson`, `ujson` to speed up serialization and deserialization.
  warnings.warn(

ts.plot(kind="box")

<Axes: >

ts.plot(kind="hist")

<Axes: ylabel='Frequency'>

# how many values from each year?
ts.index.strftime("%Y").value_counts().sort_index().plot(kind="bar")

<Axes: >

Missing Values

ts_missing = ts.copy().to_numpy()

# add missing values for demonstration
ts_missing[::10] = np.nan
ts_missing[::11] = np.nan

plt.plot(ts_missing)
plt.show()

Removing missing values

ts_cleaned = pd.Series(ts_missing).dropna()

ts_cleaned.plot()

<Axes: >

Imputing Missing Values

from sktime.transformations.series.impute import Imputer

imputer = Imputer(method="mean")

ts_imputed = imputer.fit_transform(ts_missing)

plt.plot(ts_imputed, color="C1", label="Imputed", linestyle="--")
plt.plot(ts_missing, color="C0", label="Missing")
plt.legend()
plt.show()

imputer = Imputer(method="ffill")
ts_imputed = imputer.fit_transform(ts_missing)
plt.plot(ts_imputed, color="C1", label="Imputed", linestyle="--")
plt.plot(ts_missing, color="C0", label="Missing")
plt.legend()
plt.show()

imputer = Imputer(method="bfill")
ts_imputed = imputer.fit_transform(ts_missing)
plt.plot(ts_imputed, color="C1", label="Imputed", linestyle="--")
plt.plot(ts_missing, color="C0", label="Missing")
plt.legend()
plt.show()

imputer = Imputer(method="linear")
ts_imputed = imputer.fit_transform(ts_missing)
plt.plot(ts_imputed, color="C1", label="Imputed", linestyle="--")
plt.plot(ts_missing, color="C0", label="Missing")
plt.legend()
plt.show()

Anomalies

from sktime.transformations.series.func_transform import FunctionTransformer

Hampel Filter

from sktime.transformations.series.outlier_detection import HampelFilter

# if we lower the threshold we would detect more outliers (default n_sigma is 3)
annotator = HampelFilter(window_length=12, n_sigma=3, return_bool=True)

labels = annotator.fit_transform(ts).ravel()

/var/folders/kj/v66zvn217x31k6lx63lt02q40000gn/T/ipykernel_42309/3923777818.py:1: FutureWarning: Series.ravel is deprecated. The underlying array is already 1D, so ravel is not necessary.  Use `to_numpy()` for conversion to a numpy array instead.
  labels = annotator.fit_transform(ts).ravel()

plt.plot(ts.ravel())
plt.scatter(np.arange(len(labels)), ts, alpha=labels, color="red", marker="x")
plt.show()

/var/folders/kj/v66zvn217x31k6lx63lt02q40000gn/T/ipykernel_42309/4027205443.py:1: FutureWarning: Series.ravel is deprecated. The underlying array is already 1D, so ravel is not necessary.  Use `to_numpy()` for conversion to a numpy array instead.
  plt.plot(ts.ravel())

Grubbs Test

from scipy import stats

Estimation

Each point is judged using statistics computed from the full series

def grubbs_test(ts, alpha=0.05):
    z_scores = np.abs(stats.zscore(ts))
    threshold = stats.t.ppf((1 + (1-alpha))/2, len(ts) - 2)
    return (z_scores > threshold).astype(int)

# you can use the function directly on the data
labels = grubbs_test(ts)

# or you can use the FunctionTransformer (the advantage is that this is compatible with all sktime features)
transformer = FunctionTransformer(func=grubbs_test, X_type="pd.Series")
labels = transformer.fit_transform(ts)

plt.plot(ts.to_numpy())
plt.scatter(np.arange(len(labels)), ts.to_numpy(), alpha=labels.to_numpy(), color="red", marker="x")
plt.show()

Prediction

Processes the time series causally, using only past data available up to each time point

def grubbs_test_prediction(ts, alpha=0.05, starting_win=12):
    initial_labels = list(grubbs_test(ts[:starting_win], alpha))
    for i in range(starting_win, len(ts)):
        initial_labels.append(grubbs_test(np.array(ts)[:i+1], alpha)[-1])  # -1 to only get the label of the last point
    return np.array(initial_labels)

# you can use the function directly on the data
labels = grubbs_test_prediction(ts)

# or you can use the FunctionTransformer (the advantage is that this is compatible with all sktime features)
transformer = FunctionTransformer(func=grubbs_test_prediction, X_type="pd.Series")
labels = transformer.fit_transform(ts)

plt.plot(ts.to_numpy())
plt.scatter(np.arange(len(labels)), ts.to_numpy(), alpha=labels.to_numpy(), color="red", marker="x")
plt.show()

IQR-based

Estimation

def iqr_outliers(x, k=1.5):
    q1 = np.quantile(x, 0.25)
    q3 = np.quantile(x, 0.75)
    iqr = q3 - q1
    outliers = list()
    for v in x:
        if v < q1 - k * iqr or v > q3 + k * iqr:
            outliers.append(True)
        else:
            outliers.append(False)
    return np.array(outliers).astype(int)

# again, we decrease the threshold to detect more outliers (default is 1.5)
transformer = FunctionTransformer(func=iqr_outliers, kw_args=dict(k=1), X_type="pd.Series")

labels = transformer.fit_transform(ts)
labels

0      0
1      0
2      0
3      0
4      0
      ..
139    1
140    0
141    0
142    0
143    0
Length: 144, dtype: int64

plt.plot(ts.to_numpy())
plt.scatter(np.arange(len(labels)), ts.to_numpy(), alpha=labels, color="red", marker="x")
plt.show()

Prediction

def iqr_outliers_prediction(x, k=1.5, starting_win=12):
    initial_labels = list(iqr_outliers(x[:starting_win], k))
    for i in range(starting_win, len(x)):
        initial_labels.append(iqr_outliers(np.array(x)[:i+1], k)[-1])
    return np.array(initial_labels)

# you can use the function directly on the data
labels = iqr_outliers_prediction(ts)

# or you can use the FunctionTransformer (the advantage is that this is compatible with all sktime features)
transformer = FunctionTransformer(func=iqr_outliers_prediction, kw_args=dict(k=1.5), X_type="pd.Series")
labels = transformer.fit_transform(ts)

plt.plot(ts.to_numpy())
plt.scatter(np.arange(len(labels)), ts.to_numpy(), alpha=labels, color="red", marker="x")
plt.show()

Replacing Anomalies

transformer = FunctionTransformer(func=iqr_outliers, kw_args=dict(k=1), X_type="pd.Series")
labels = transformer.fit_transform(ts)

plt.plot(ts.to_numpy())
plt.scatter(np.arange(len(labels)), ts.to_numpy(), alpha=labels, color="red", marker="x")
plt.show()

ts_replaced = ts.copy().to_numpy()
ts_replaced[labels.astype(bool)] = np.nan

plt.plot(ts_replaced)
plt.show()

imputer = Imputer(method="nearest")

ts_imputed = imputer.fit_transform(ts_replaced)

plt.plot(ts.to_numpy(), label="Original")
plt.plot(ts_imputed, label="Imputed")
plt.legend()
plt.show()

Manipulating Time Series

# let's split the time series into two
ts1 = ts[0: len(ts)//2].to_numpy()
ts2 = ts[len(ts)//2:].to_numpy()

plt.plot(ts1)
plt.plot(ts2)
plt.show()

plt.plot(ts1)
plt.plot(ts2)
plt.yscale("log")
plt.show()

Normalizations

from sktime.transformations.series.adapt import TabularToSeriesAdaptor

Min-Max normalization

from sklearn.preprocessing import MinMaxScaler

scaler = TabularToSeriesAdaptor(MinMaxScaler(), fit_in_transform=True)

ts1_minmax = scaler.fit_transform(ts1.reshape(-1, 1))
ts2_minmax = scaler.fit_transform(ts2.reshape(-1, 1))

plt.plot(ts1_minmax)
plt.plot(ts2_minmax)
plt.axhline(0, color="black", linestyle="--")
plt.axhline(1, color="black", linestyle="--")
plt.show()

Offset Transformation

from sklearn.preprocessing import StandardScaler

scaler = TabularToSeriesAdaptor(StandardScaler(with_std=False), fit_in_transform=True)

ts1_offset = scaler.fit_transform(ts1.reshape(-1, 1))
ts2_offset = scaler.fit_transform(ts2.reshape(-1, 1))

ts1_offset.mean(), ts1_offset.std()

(6.315935428978668e-15, 47.371803533977136)

ts2_offset.mean(), ts2_offset.std()

(-1.2631870857957336e-14, 85.8368347002499)

plt.plot(ts1_offset)
plt.plot(ts2_offset)
plt.show()

Z-score normalization (Amplitude Scaling)

scaler = TabularToSeriesAdaptor(StandardScaler(), fit_in_transform=True)

ts1_zscored = scaler.fit_transform(ts1.reshape(-1, 1))
ts2_zscored = scaler.fit_transform(ts2.reshape(-1, 1))

ts1_zscored.mean(), ts1_zscored.std()

(1.4494578377050655e-16, 0.9999999999999999)

ts2_zscored.mean(), ts2_zscored.std()

(-1.603655480014115e-16, 1.0)

plt.plot(ts1_zscored)
plt.plot(ts2_zscored)
plt.show()

Smoothing

def moving_average(x, w):
    return x.rolling(window=w).mean()

transformer = FunctionTransformer(func=moving_average, kw_args=dict(w=4), X_type="pd.Series")

ts1_smoothed = transformer.fit_transform(pd.Series(ts1_zscored.ravel()))
ts2_smoothed = transformer.fit_transform(pd.Series(ts2_zscored.ravel()))

plt.plot(ts1_zscored.ravel())
plt.plot(ts1_smoothed)

[<matplotlib.lines.Line2D at 0x179fdb560>]

plt.plot(ts2_zscored.ravel())
plt.plot(ts2_smoothed)

[<matplotlib.lines.Line2D at 0x179ddce90>]

plt.plot(ts1_smoothed)
plt.plot(ts2_smoothed)
plt.show()

Log transformation

from sktime.transformations.series.boxcox import LogTransformer

log_transformer = LogTransformer()

ts1_log = log_transformer.fit_transform(ts1)
ts2_log = log_transformer.fit_transform(ts2)

plt.plot(ts1_log)
plt.plot(ts2_log)

[<matplotlib.lines.Line2D at 0x17a023a40>]

Differencing

from sktime.transformations.series.difference import Differencer

differencer =  Differencer(1)

ts1_diff = differencer.fit_transform(ts1)
ts2_diff = differencer.fit_transform(ts2)

plt.plot(ts1_diff)
plt.plot(ts2_diff)
plt.show()

Components

from sktime.forecasting.trend import STLForecaster

forecaster = STLForecaster(sp=12)

forecaster.fit(ts1)
ts1_trend = forecaster.trend_
ts1_seasonal = forecaster.seasonal_
ts1_resid = forecaster.resid_

forecaster.fit(ts2)
ts2_trend = forecaster.trend_
ts2_seasonal = forecaster.seasonal_
ts2_resid = forecaster.resid_

fig, axs = plt.subplots(4, 1, sharex=True)
axs[0].plot(ts1, color="C0")
axs[0].set_ylabel("Original")
axs[1].plot(ts1_trend, color="C1")
axs[1].set_ylabel("Trend")
axs[2].plot(ts1_seasonal, color="C2")
axs[2].set_ylabel("Seasonal")
axs[3].plot(ts1_resid, color="C3")
axs[3].set_ylabel("Residual")
plt.show()

fig, axs = plt.subplots(4, 1,  sharex=True)
axs[0].plot(ts2, color="C0")
axs[0].set_ylabel("Original")
axs[1].plot(ts2_trend, color="C1")
axs[1].set_ylabel("Trend")
axs[2].plot(ts2_seasonal, color="C2")
axs[2].set_ylabel("Seasonal")
axs[3].plot(ts2_resid, color="C3")
axs[3].set_ylabel("Residual")
plt.show()

Stationarity

from sktime.param_est.stationarity import StationarityADF
from statsmodels.tsa.stattools import adfuller

adfuller(ts)

(0.8153688792060371,
 0.9918802434376408,
 13,
 130,
 {'1%': -3.4816817173418295,
  '5%': -2.8840418343195267,
  '10%': -2.578770059171598},
 996.692930839019)

sty_est = StationarityADF()  
sty_est.fit(ts)  
sty_est.stationary_

False

Detrend

from sktime.transformations.series.detrend import Detrender

# linear detrending
detrender = Detrender()

ts_detrended = detrender.fit_transform(ts)

plt.plot(ts_detrended.to_numpy())
plt.show()

Deseasonalize

from sktime.transformations.series.detrend import Deseasonalizer
from sktime.utils.plotting import plot_correlations

plot_correlations(ts)
plt.show()

deseasonalizer = Deseasonalizer(sp=12, model="multiplicative")  # yearly seasonality

ts_deseasonalized = deseasonalizer.fit_transform(ts)

plt.plot(ts_deseasonalized.to_numpy())

[<matplotlib.lines.Line2D at 0x16c016c00>]

Detrend and Deseasonalize

detr_dese = Detrender(model="multiplicative") * Deseasonalizer(sp=12, model="multiplicative")  # this is a pipeline
detr_dese

TransformerPipeline(steps=[Detrender(model='multiplicative'),
                           Deseasonalizer(model='multiplicative', sp=12)])

Please rerun this cell to show the HTML repr or trust the notebook.

ts_detr_dese = detr_dese.fit_transform(ts)

sty_est = StationarityADF()  
sty_est.fit(ts_detr_dese)  
sty_est.stationary_

True

plt.plot(ts_detr_dese.to_numpy())

[<matplotlib.lines.Line2D at 0x1797e5b80>]

Pipeline

# without pipeline
ts1_log = LogTransformer().fit_transform(ts1)
ts1_det = Detrender().fit_transform(ts1_log)
ts1_scaled = StandardScaler().fit_transform(ts1_det)

from sktime.pipeline import make_pipeline

pipe = make_pipeline(LogTransformer(), Detrender(), TabularToSeriesAdaptor(StandardScaler(), fit_in_transform=True))

ts1_pipe = pipe.fit_transform(ts1.reshape(-1, 1))
ts2_pipe = pipe.fit_transform(ts2.reshape(-1, 1))

plt.plot(ts1_pipe)
plt.plot(ts2_pipe)
plt.show()

Multivariate Time Series

from sktime.datasets import load_uschange

_, X = load_uschange()

sns.lineplot(X[["Savings", "Income"]])
plt.show()

pipe = make_pipeline(TabularToSeriesAdaptor(MinMaxScaler(), fit_in_transform=True), FunctionTransformer(func=np.log1p), Differencer(1))
# this is equivalent to the following
# pipe = TabularToSeriesAdaptor(MinMaxScaler(), fit_in_transform=True) * FunctionTransformer(func=np.log1p) * Differencer(1)
pipe

TransformerPipeline(steps=[TabularToSeriesAdaptor(fit_in_transform=True,
                                                  transformer=MinMaxScaler()),
                           FunctionTransformer(func=<ufunc 'log1p'>),
                           Differencer()])

Please rerun this cell to show the HTML repr or trust the notebook.

X_transformed = pipe.fit_transform(X)

sns.lineplot(X_transformed[["Savings", "Income"]])

<Axes: >