<jupyter_start><jupyter_text>ARIMAAn [AutoRegressive Integrated Moving Average](https://en.wikipedia.org/wiki/Autoregressive_integrated_moving_average) model is a popular model used in time series analysis to understand the data or forecast future points.This implementation can fit a model to each time series in a batch and perform in-sample predictions and out-of-sample forecasts. It is designed to give the best performance when working on a large batch of time series.Useful links:- cuDF documentation: https://docs.rapids.ai/api/cudf/stable- cuML's ARIMA API docs: https://docs.rapids.ai/api/cuml/stable/api.htmlarima- a good introduction to ARIMA: https://otexts.com/fpp2/arima.html  Setup Imports<jupyter_code>import cudf
from cuml.tsa.arima import ARIMA

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import random<jupyter_output><empty_output><jupyter_text>Loading utilThe data for this demo is stored in a simple CSV format:- the data series are stored in columns- the first column contains the date of the data points- the first row contains the name of each variableFor example, let's check the *population estimate* dataset:<jupyter_code>!cat data/time_series/population_estimate.csv | head<jupyter_output><empty_output><jupyter_text>We define a helper function to load a dataset with a given name and return a GPU dataframe. We discard the date, and limit the batch size for convenience.<jupyter_code>def load_dataset(name, max_batch=4):
    import os
    pdf = pd.read_csv(os.path.join("data", "time_series", "%s.csv" % name))
    return cudf.from_pandas(pdf[pdf.columns[1:max_batch+1]].astype(np.float64))<jupyter_output><empty_output><jupyter_text>Visualization utilWe define a helper function that displays the data, and optionally a prediction starting from a given index. Each time series is plot separately for better readability.<jupyter_code>def visualize(y, pred=None, pred_start=None, lower=None, upper=None):
    n_obs, batch_size = y.shape
    col = ["#1f77b4", "#ff7f0e"]

    # Create the subplots
    c = min(batch_size, 2)
    r = (batch_size + c - 1) // c
    fig, ax = plt.subplots(r, c, squeeze=False)
    ax = ax.flatten()
    
    # Range for the prediction
    if pred is not None:
        pred_start = n_obs if pred_start is None else pred_start
        pred_end = pred_start + pred.shape[0]
    else:
        pred_end = n_obs
    
    # Plot the data
    for i in range(batch_size):
        title = y.columns[i]
        if pred is not None:
            ax[i].plot(np.r_[pred_start:pred_end],
                       pred[pred.columns[i]].to_numpy(),
                       linestyle="--", color=col[1])
        # Prediction intervals
        if lower is not None and upper is not None:
            ax[i].fill_between(np.r_[pred_start:pred_end],
                               lower[lower.columns[i]].to_numpy(),
                               upper[upper.columns[i]].to_numpy(),
                               alpha=0.2, color=col[1])
        ax[i].plot(np.r_[:n_obs], y[title].to_numpy(), color=col[0])
        ax[i].title.set_text(title)
        ax[i].set_xlim((0, pred_end))
    for i in range(batch_size, r*c):
        fig.delaxes(ax[i])
    fig.tight_layout()
    fig.patch.set_facecolor('white')
    plt.show()<jupyter_output><empty_output><jupyter_text>Non-seasonal ARIMA modelsA basic `ARIMA(p,d,q)` model is made of three components: - An **Integrated** (I) component: the series is differenced `d` times until it is stationary - An **AutoRegressive** (AR) component: the variable is regressed on its `p` past values - A **Moving Average** (MA) component: the variable is regressed on `q` past error termsThe model can also incorporate an optional constant term (called *intercept*). A simple MA(2) exampleWe start with a simple Moving Average model. Let's first load and visualize the *migrations in Auckland by age* dataset:<jupyter_code>df_mig = load_dataset("net_migrations_auckland_by_age", 4)
visualize(df_mig)<jupyter_output><empty_output><jupyter_text>We want to fit the model with `q`=2 and with an intercept.The `ARIMA` class accepts cuDF dataframes or array-like types as input (host or device), e.g numpy arrays. Here we already have a dataframe so we can simply pass it to the `ARIMA` constructor with the model parameters:<jupyter_code>model_mig = ARIMA(df_mig, order=(0,0,2), fit_intercept=True)
model_mig.fit()<jupyter_output><empty_output><jupyter_text>We can now forecast and visualize the results:<jupyter_code>fc_mig = model_mig.forecast(10)
visualize(df_mig, fc_mig)<jupyter_output><empty_output><jupyter_text>If we want to get the parameters that were fitted to the model, we can use `get_fit_params` or the corresponding properties. The parameters are organized in 2D arrays: one row represents one parameter and the columns are different batch members.<jupyter_code>print(model_mig.get_fit_params()["ma"])
print(model_mig.ma_)<jupyter_output><empty_output><jupyter_text>It is also possible to get a compact numpy array containing all the parameters with `pack`, or similarly to load the parameters into a model with `unpack`:<jupyter_code>print(model_mig.pack())<jupyter_output><empty_output><jupyter_text>We can also get the log-likelihood of the parameters w.r.t to the series, and evaluate various information criteria:<jupyter_code>print("log-likelihood:\n", model_mig.llf)
print("\nAkaike Information Criterion (AIC):\n", model_mig.aic)
print("\nCorrected Akaike Information Criterion (AICc):\n", model_mig.aicc)
print("\nBayesian Information Criterion (BIC):\n", model_mig.bic)<jupyter_output><empty_output><jupyter_text>An ARIMA(1,2,1) exampleLet's now load the *population estimate* dataset. For this dataset a first difference is not enough to make the data stationary because of the quadratic trend, so we decide to go with `d`=2.This time we won't simply forecast but also predict in-sample:<jupyter_code>df_pop = load_dataset("population_estimate")

# Fit an ARIMA(1,2,1) model
model_pop = ARIMA(df_pop, order=(1,2,1), fit_intercept=True)
model_pop.fit()

# Predict in-sample and forecast out-of-sample
pred_pop = model_pop.predict(80, 160)
visualize(df_pop, pred_pop, 80)<jupyter_output><empty_output><jupyter_text>Confidence intervalsTo get confidence intervals when forecasting, we can specify the confidence level (here 95%):<jupyter_code>fc_pop, lower_pop, upper_pop = model_pop.forecast(23, level=0.95)
visualize(df_pop, fc_pop, lower=lower_pop, upper=upper_pop)<jupyter_output><empty_output><jupyter_text>Seasonal ARIMA models[Seasonal ARIMA models](https://otexts.com/fpp2/seasonal-arima.html) are expressed in the form `ARIMA(p,d,q)(P,D,Q)s` and have additional seasonal components that we denote SAR and SMA.We can also choose to apply a first or second seasonal difference, or combine a non-seasonal and a seasonal difference (note: `p+P <= 2` is required). An ARIMA(1,1,1)(1,1,1)12 exampleWe load the *guest nights by region* dataset. This dataset shows a strong seasonal component with a period of 12 (annual cycle, monthly data), and also a non-seasonal trend. A good choice is to go with `d`=1, `D`=1 and `s`=12.We create the model with seasonal parameters, and forecast:<jupyter_code>df_guests = load_dataset("guest_nights_by_region", 4)

# Create and fit an ARIMA(1,1,1)(1,1,1)12 model:
model_guests = ARIMA(df_guests, order=(1,1,1), seasonal_order=(1,1,1,12),
                     fit_intercept=False)
model_guests.fit()
# Forecast
fc_guests = model_guests.forecast(40)

# Visualize after the time step 200
visualize(df_guests[200:], fc_guests)<jupyter_output><empty_output><jupyter_text>Missing observationsARIMA supports missing observations in the data. You can also pad your dataset at the start in order to batch computations even if the series have different lengths.To illustrate that, let's create a fake dataset from the seasonal dataset used above. We will simulate series of different lengths and add missing observations:<jupyter_code># Cut dataset to 100 observations
df_guests_missing = df_guests[:100].copy()

for title in df_guests_missing.columns:
    # Missing observations at the start to simulate varying lengths
    n_leading = random.randint(5, 40)
    df_guests_missing[title][:n_leading]=None
    
    # Random missing observations in the middle
    missing_obs = random.choices(range(n_leading, 100), k=random.randint(5, 20))
    df_guests_missing[title][missing_obs]=None<jupyter_output><empty_output><jupyter_text>Note that missing observations need to be represented by the value `NaN` to convert the dataset to a numeric array. `NA`s in dataframes can be filled with:<jupyter_code>df_guests_missing = df_guests_missing.fillna(np.nan)<jupyter_output><empty_output><jupyter_text>We can now fit a model. Here we will do in- and out-of-sample predictions, to showcase how this model can fill the gaps in data.<jupyter_code># Create and fit an ARIMA(1,1,1)(1,1,1)12 model:
model_guests_missing = ARIMA(df_guests_missing, order=(1,1,1), seasonal_order=(1,1,1,12),
                             fit_intercept=False)
model_guests_missing.fit()
# Forecast
fc_guests_missing = model_guests_missing.predict(0, 120)

visualize(df_guests_missing, fc_guests_missing, 0)<jupyter_output><empty_output><jupyter_text>Note that the model can't form predictions at the start where we padded with missing values. The first in-sample predictions will be equal to a constant value (0 in the absence of intercept).  Exogenous variablesARIMA also supports exogenous variables. As this model works with a batch, there are some limitations: each column of `endog` corresponds to a fixed number of columns of `exog`. Exogenous variables can't be shared by different series (they have to be duplicated if necessary), and all series must have the same number of exogenous variables. The shape of `exog` must be `(n_obs, batch_size * n_exog)`, and columns are grouped by corresponding batch id. For predictions and forecasts, values of the exogenous variables for future steps must be provided in an array of shape `(nsteps, batch_size * n_exog)`.Note that endogenous variables might contain missing observations but exogenous variables cannot.To illustrate this, we will again create a fake dataset from the one used above, adding some procedural variables.<jupyter_code>nb = 4

# Generate exogenous variables and coefficients
get_sine = lambda n, period: \
    np.sin(np.r_[:n] * 2 * np.pi / period + np.random.uniform(0, period))
np_exog = np.column_stack([get_sine(319, T)
                           for T in np.random.uniform(20, 100, 2 * nb)])
np_exog_coef = np.random.uniform(20, 200, 2 * nb)

# Create dataframes for the past and future values
df_exog = cudf.DataFrame(np_exog[:279])
df_exog_fut = cudf.DataFrame(np_exog[279:])

# Add linear combination of the exogenous variables to the endogenous
df_guests_exog = df_guests.copy()
for ib in range(nb):
    df_guests_exog[df_guests_exog.columns[ib]] += \
        np.matmul(np_exog[:279, ib*2:(ib+1)*2], np_exog_coef[ib*2:(ib+1)*2])
# Create and fit an ARIMA(1,0,1)(1,1,1)12 (c) model with exogenous variables
model_guests_exog = ARIMA(endog=df_guests_exog, exog=df_exog,
                          order=(1,0,1), seasonal_order=(1,1,1,12),
                          fit_intercept=True)
model_guests_exog.fit()
# Forecast
fc_guests_exog = model_guests_exog.forecast(40, exog=df_exog_fut)

# Visualize after the time step 100
visualize(df_guests_exog[100:], fc_guests_exog)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09474050853361198, 'fraction_numerical': 0.016457680250783698, 'mean_word_length': 4.838840874428063, 'pattern_counts': {'":': 0, '<': 70, '<?xml version=': 0, '>': 69, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '370616', 'license': 'Apache-2.0', 'n_tokens_mistral': 3833, 'n_tokens_neox': 3642, 'n_words': 1171, 'path': 'notebooks/arima_demo.ipynb', 'repo_name': 'daxiongshu/cuml', 'size': '11484'}