Time series analysis with Granite Time

!pip install "granite-tsfm[notebooks] @ git+https://github.com/ibm-granite/granite-tsfm.git@v0.3.1";

Next, import the libraries that contain and support TSPulse:

from tsfm_public.models.tspulse import TSPulseForReconstruction
from tsfm_public.models.tspulse.utils.helpers import get_embeddings
from tsfm_public.models.tspulse import TSPulseForClassification
from tsfm_public.toolkit.dataset import ClassificationDFDataset
from tsfm_public import TimeSeriesPreprocessor
from tsfm_public.models.tspulse import TSPulseForReconstruction
from tsfm_public.toolkit.time_series_imputation_pipeline import TimeSeriesImputationPipeline
from tsfm_public.toolkit.time_series_classification_preprocessor import (
    TimeSeriesClassificationPreprocessor,
)
from tsfm_public.models.tspulse.modeling_tspulse import TSPulseForReconstruction
from tsfm_public.toolkit.ad_helpers import AnomalyScoreMethods
from tsfm_public.toolkit.time_series_anomaly_detection_pipeline import TimeSeriesAnomalyDetectionPipeline
from tsfm_public.toolkit.util import convert_tsfile_to_dataframe
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import utils
import torch

The first time series dataset is an example of time series sales data that contains the sales for each day across 3 different shops.

df = pd.read_csv("ice_cream_sales.csv", index_col=False)
df.head(5)

This should output:

Unnamed: 0	date	Pike Place Sales	Fremont Sales	Ballard Sales	Event	Festival	Temp Max	Temp Min	Precip
0	1	2022-01-01	1553.173214	1144.518000	1758.323182	True	True	34	16	0.00
1	2	2022-01-02	1211.258185	739.110199	1244.401844	False	False	38	28	1.30
2	3	2022-01-03	1342.835436	820.118584	1342.835436	False	False	42	33	0.52
3	4	2022-01-04	1442.806029	1065.528522	1842.806029	False	True	42	34	0.05
4	5	2022-01-05	1431.781942	875.601802	1431.781942	False	False	39	35	0.06
5	6	2022-01-06	1291.068914	776.108094	1291.068914	False	False	51	36	1.06

The data has an extra column that you’ll drop by using the drop() method. Next, set the date column to a Pandas DateTime object:

df = df.drop("Unnamed: 0", axis=1)
df_with_nans = df
df_with_nans['date'] = pd.to_datetime(df_with_nans['date'])

One feature of TSPulse is that, unlike more statistical algorithms, it performs well on nonstationary time series data. This means that you don’t have to check for stationarity, decompose the time series or perform differencing. Linear regression models usually require that you check the autocorrelation function (ACF) and partial autocorrelation function (PACF) to determine which lags are most predictive. With TSPulse, we can often simply skip preprocessing and work with multiple time series datasets immediately.

timestamp_column = "date"
id_columns = []  # mention the ids that uniquely identify a time-series.
target_columns = [
    "Pike Place Sales",
    "Fremont Sales",
    "Ballard Sales"
]  # mention the target column names in the dataset that should be imputed by the model

print(df_with_nans.shape)
print(df_with_nans.head())

column_specifiers = {
    "timestamp_column": timestamp_column,
    "id_columns": id_columns,
    "target_columns": target_columns,
    "control_columns": [] # does adding these help?
}

Now, create a visualization of the fluctuations in sales data:

plt.figure(figsize=(16,6))
plt.plot(df_with_nans['date'], df_with_nans["Pike Place Sales"])
plt.title("Pike Place Sales")

We can see that the time series has a strong seasonal component (often called “seasonality”) and a weak stochastic variation across days.

To see how TSPulse can help the owner of the ice cream shops, imagine that they lose 10% of their sales data. The missing values create nonlinear relationships from one period of time to another, making forecasting tricky. To create those missing values a utility function, which is available in the GitHub repository, will set about 10% of the sales values to NaN. About half of the NaN values will be individual data points and about half will be a sequence of data points.

import random

df_with_nans = utils.introduce_nans(df_with_nans, targets=target_columns, nan_fraction=0.1)
df_with_nans.head(10)

The imputation tracks the trends in a multivariate time series and uses them to extrapolate what missing values might be. This task is common in time series analysis because data gets corrupted, errors happen in data collection and pipelines break down.

Now create a TimeSeriesPreprocessor to the target columns in the data.

tsp = TimeSeriesPreprocessor(
    **column_specifiers,
    context_length=512,
    prediction_length=0,
    scaling=True,
    encode_categorical=False,
    scaler_type="standard",
)

tsp.train(df_with_nans)  # train the tsp

Now download the TSPulse model itself and create a pipeline to use it for imputation:

model = TSPulseForReconstruction.from_pretrained(
    "ibm-granite/granite-timeseries-tspulse-r1",
    revision="tspulse-hybrid-dualhead-512-p8-r1",
    num_input_channels=tsp.num_input_channels,
    mask_type="user", 
)

device = "cuda" if torch.cuda.is_available() else "cpu"

pipe = TimeSeriesImputationPipeline(model, feature_extractor=tsp, batch_size=1000, device=device)
imputed_df = pipe(df_with_nans)

Graph the output of the model against the ground truth:

plt.figure(figsize=(15, 9))

number_of_days_to_plot = 365

#x_range = np.arange(0, len(imputed_df), dtype=float)
x_range = np.arange(0, number_of_days_to_plot, dtype=float)

for i, base_col in enumerate(target_columns):

    imputed_col = f"{base_col}_imputed"
    print(imputed_col)
    observed_vals = df_with_nans[base_col][:number_of_days_to_plot]
    imputed_vals = imputed_df[imputed_col][:number_of_days_to_plot]

    pos_observed = ~observed_vals.isna()
    pos_imputed = observed_vals.isna()
    plt.subplot(len(target_columns), 1, i + 1)

    full_vals = df[base_col][:number_of_days_to_plot]
    y_min = np.min(full_vals)
    y_max = np.max(full_vals)

    plt.vlines(x=x_range[pos_imputed], ymin=y_min, ymax=y_max, color='#f2f2f2', linestyle='-', zorder=0)
    plt.plot(x_range, full_vals, color="blue", linewidth=1, label="Ground_Truth")  # actual ground truth
    plt.scatter(x_range[pos_imputed], imputed_vals[pos_imputed], color="red", s=5, label="Imputed")

    plt.title(base_col)
    plt.legend()
    plt.grid(False)

plt.tight_layout()
plt.show()

To get a precise accuracy metric, you’ll calculate the RMSE, MAE and MAPE metrics

actual = df[target_columns].to_numpy(dtype=float)
imputed_values = imputed_df[[f"{col}_imputed" for col in target_columns]]  # df having imputed values at the missing data positions
imputed = imputed_values.to_numpy(dtype=float)

missing_positions = np.isnan(df_with_nans[target_columns])

print(pd.DataFrame(
    {
        "results": {
            "root_mean_squared_error": np.sqrt(np.mean(np.square(actual[missing_positions] - imputed[missing_positions]))),
            "mean_absolute_error": np.mean(np.abs(actual[missing_positions] - imputed[missing_positions])),
            "mape": np.mean(np.abs((actual[missing_positions] - imputed[missing_positions]) / actual[missing_positions])) * 100,
        }
    }
))

This will output:

                            results
root_mean_squared_error  106.868017
mean_absolute_error       75.099967
mape                       5.273472

The imputed values are only 5% off from the actual values. Not perfect, but enough to save the shop from struggling with lost data.

The next task is to look at anomalies and outliers that diverge from the usual underlying patterns and seasonal variations. In many domains, an anomaly can be an early warning of trouble. For instance, a sensor spike in industrial machinery might indicate failing equipment or a sudden drop in website traffic might indicate server issues. However, not all anomalies are bad. Sometimes anomalies are interesting and worth accounting for in planning and forecasting. A viral marketing campaign can cause unusual sales growth or a special event might cause a traffic surge on city streets. In either case, having a tool to detect anomalous signals in a time series is a powerful tool. The ice cream shops are interested in what days had anomalous sales so that they can explore what might have caused those anomalies. To do this, you’ll examine one location by isolating the “Ballard Sales” column. To do this, create a TimeSeriesAnomalyDetectionPipeline and set the target_colums to “Ballard Sales”.

anomaly_model = TSPulseForReconstruction.from_pretrained(
    "ibm-granite/granite-timeseries-tspulse-r1",
    num_input_channels=1,
    revision="main",
    mask_type="user",
)

pipeline = TimeSeriesAnomalyDetectionPipeline(
    anomaly_model,
    timestamp_column="date",
    target_columns=["Ballard Sales"],
    prediction_mode=[
        AnomalyScoreMethods.FREQUENCY_RECONSTRUCTION.value,
        AnomalyScoreMethods.PREDICTIVE.value
    ],
    aggregation_length=32,
    aggr_function="max",
    smoothing_length=1, #no smoothing, we're looking for 1 off events
    least_significant_scale=0.01,
    least_significant_score=0.1,
)

result = pipeline(df, batch_size=32, predictive_score_smoothing=False)

The anomaly detection needs 512 data points of past values to set a baseline against which to detect anomalies, so visualize the results of the anomaly detection after the 512th day:

fig, ax = plt.subplots(1, 1, figsize=(15, 3))
ax2 = ax.twinx()

date_array  = df[512:]['date']

bsales = ax.plot(date_array, result[512:]['Ballard Sales'], linewidth=1, linestyle = 'solid', label="Ballard Sales")
anom = ax2.plot(date_array, result[512:].anomaly_score, linewidth=1, color="orange", label="Anomaly Score")

ax.legend(loc=2)
ax2.legend(loc=1)

ax.set_title("Anomaly Detection after 512th date", fontsize=16)

The anomalous values are mostly on relatively warm and dry winter days:

result[result["anomaly_score"] > 0.9]

This will output:

date	Pike Place Sales	Fremont Sales	Ballard Sales	Event	Festival	Temp Max	Temp Min	Precip	anomaly_score
609	2023-09-02	1796.424619	1244.250202	2209.256063	False	True	84	58	0.00	0.994726
624	2023-09-17	1744.949799	1062.452261	1750.369403	False	False	70	55	0.00	0.924997
679	2023-11-11	1821.270839	1317.356317	2036.171961	True	True	54	46	0.20	0.991932
1043	2024-11-09	1644.689543	1208.192755	2065.722624	False	True	53	48	0.02	0.991474

Finally, TSPulse offers functionality to classify time series data as well. This technique is a powerful way to analyze a sequence of data and classify the sequence. The ice cream shops have kept track of the number of customers in the store during special events so they can make informed decisions about staffing and scheduling. The dataset consists of 200 events with type of event and number of customers in the store averaged across a 10-minute increment. With this exercise, the ice cream shop wants to determine the pattern of customer visits based on the type of event happening around the store.

sales_patterns = pd.read_csv("ic_patterns.csv")
sales_patterns = sales_patterns.drop("Unnamed: 0", axis=1)
sales_patterns.head(5)

The event data needs to be formatted so that each row contains the class of the event with a Pandas Series that contains all the readings for the time series. To create this, take a list of events and create a row with the type of the event and all the measurements for that event.

events = []
event_ids = sales_patterns["event"].unique()
for e in event_ids:
    events.append({"type":sales_patterns[sales_patterns["event"] == e].iloc[0]["type"], "ts":pd.Series(sales_patterns[sales_patterns["event"] == e]["sales"].tolist())})

event_time_series = pd.DataFrame(events)

# plot the first event:
event_time_series.iloc[0]['ts'].plot()

Now you’re ready to create two datasets: a test set and a train set.

test_dataset, train_dataset = utils.test_train_dataset(event_time_series, test_size = 0.25, input_columns = ['ts'], label_column=['type'])

With the datasets ready, download the classifier from HuggingFace and configure it for our dataset and use case.

config_dict = {
    "head_reduce_d_model": 1,
    "decoder_mode": "mix_channel",
    "head_gated_attention_activation": "softmax",
    "mask_ratio": 0.3,
    "channel_virtual_expand_scale": 1,
    "loss": "cross_entropy",
    "disable_mask_in_classification_eval": True,
    "ignore_mismatched_sizes": True,
    "num_input_channels" : 1, # how many different series are in each time series
    "num_targets" : 4 # how many classes are in the dataset?
}

classifier_model = TSPulseForClassification.from_pretrained(
    "ibm-granite/granite-timeseries-tspulse-r1", revision="tspulse-block-dualhead-512-p16-r1", **config_dict
);

Before fine-tuning the model, freeze the parameters and set only the last layers to requires_grad = True. This step will update the final patch embedding layers while leaving the rest of the model untouched.

# Freezing the Backbone
for param in classifier_model.backbone.parameters():
    param.requires_grad = False

# Unfreezing the patch embedding layers
for param in classifier_model.backbone.time_encoding.parameters():
    param.requires_grad = True
for param in classifier_model.backbone.fft_encoding.parameters():
    param.requires_grad = True

The model fine-tuning process requires a validation dataset to use in training, so take 25% of the records from the train_dataset and use them for validation:

from torch.utils.data import DataLoader, random_split

valid_size = int(len(train_dataset) * 0.25)
print("validation dataset size " + str(valid_size))
train_dataset, valid_dataset = utils.random_split(train_dataset, [len(train_dataset) - valid_size, valid_size])

Now you're ready to train the model by creating a Trainer instance and passing it the model and data for fine-tuning:

finetune_trainer = utils.create_classify_trainer(classifier_model, train_dataset, valid_dataset)
finetune_trainer.train()

This process might take up to 30 minutes depending on the hardware that it’s being trained on but when it’s done, the Trainer instance contains the fine-tuned model. To test it, pass the dataset and get the most likely class for each row:

predictions_dict = finetune_trainer.predict(test_dataset)
preds_np = predictions_dict.predictions[0]

test_dataloader = DataLoader(test_dataset, batch_size=32, shuffle=False)
target_list = []
for batch in test_dataloader:
    batch_labels = batch["target_values"].numpy()
    target_list.append(batch_labels)
targets_np = np.concatenate(target_list, axis=0)

test_accuracy = np.mean(targets_np == np.argmax(preds_np, axis=1))
print("test_accuracy : ", test_accuracy * 100.0)

You can see the accuracy output:

test_accuracy :  100.0 %

While classifying real-world data might be more difficult, the accuracy that TSPulse demonstrates on benchmark datasets means that you might want to test it out for any classification task.

TSPulse is a versatile model that enables many different time series analysis tasks and can be run on limited resources like a laptop or a lightweight virtual machine.