Exercises/2_forecasting.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wb47UH9X0gqf"
      },
      "source": [
        "# University Certificate in Artificial Intelligence (Hands on AI, Third Challenge, 2022-2023, UMONS)\n",
        "# Forecasting methods\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "JflAHy-u-6EH",
        "outputId": "6799170c-8484-469f-c0da-988419adec8c"
      },
      "outputs": [],
      "source": [
        "# The dataset and package \"main\" are in the GitHub repository\n",
        "!git clone https://github.com/bsouhaib/Hands-On-AI-2022-Challenge3.git\n",
        "%cd Hands-On-AI-2022-Challenge3/Exercises\n",
        "# These packages are not installed by default\n",
        "!pip install pmdarima supersmoother"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "6MHFvCBa0gqi"
      },
      "outputs": [],
      "source": [
        "import pandas as pd\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "\n",
        "plt.rcParams[\"figure.figsize\"] = [10, 5]\n",
        "sns.set_theme()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HdVADgE70gqi"
      },
      "source": [
        "## Simulated time series (ARMA)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "msh7MtMv0gqj"
      },
      "source": [
        "We will simulate a time series from the following ARMA process:\n",
        "$$\n",
        "y_t = 0.75 y_{t-1} + 0.25 y_{t-2} + 0.65 \\varepsilon_{t-1} + 0.35 \\varepsilon_{t-2} + \\varepsilon_t\n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 339
        },
        "id": "UILM77FR0gqj",
        "outputId": "f06479bb-bcdb-42d0-f327-a990473319c1"
      },
      "outputs": [],
      "source": [
        "import statsmodels.api as sm\n",
        "\n",
        "np.random.seed(12345)\n",
        "\n",
        "T = 200\n",
        "\n",
        "# ARMA parameters\n",
        "arparams = np.array([0.75, -0.25])\n",
        "maparams = np.array([0.65, 0.35])\n",
        "ar_term = np.r_[1, -arparams]  # add zero-lag and negate\n",
        "ma_term = np.r_[1, maparams]  # add zero-lag\n",
        "\n",
        "arma_process = sm.tsa.ArmaProcess(ar_term, ma_term)\n",
        "\n",
        "y = arma_process.generate_sample(T)\n",
        "plt.plot(y)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bCSRrOrh0gqj"
      },
      "source": [
        "* Plot the autocorrelation function (ACF) of the simulated series for the first 10 lags.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 659
        },
        "id": "1mly80Kh0gqk",
        "outputId": "7efd7897-3db3-4aaf-da1d-0d96196df2a9"
      },
      "outputs": [],
      "source": [
        "from statsmodels.graphics.tsaplots import plot_acf\n",
        "\n",
        "# Hint: use plot_acf\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7gULtPy20gqk"
      },
      "source": [
        "* Fit an ARIMA model to the simulated time series using auto_arima. Do you recover the true parameters? To obtain details about your model fit, you can use the following functions: summary(), arparams(), and maparams()."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "VKB2OhD50gqk"
      },
      "outputs": [],
      "source": [
        "from pmdarima.arima import auto_arima\n",
        "\n",
        "model = ?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ADgUhlgD0gql"
      },
      "outputs": [],
      "source": [
        "def result_table(model):\n",
        "    return pd.DataFrame({\n",
        "        'Fitted parameters': [str(model.arparams()), str(model.maparams())],\n",
        "        'True parameters': [-ar_term[1:], ma_term[1:]],\n",
        "    }, index=['AR', 'MA'])\n",
        "\n",
        "result_table(model)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "1S1RmKNh0gql"
      },
      "outputs": [],
      "source": [
        "model.summary()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4WTrXT9I0gql"
      },
      "source": [
        "* Increase the size of the simulated series (e.g. $T=1000$) and refit an ARIMA model using auto_arima. Do you recover the true parameters?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "FUj4be210gql"
      },
      "outputs": [],
      "source": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NETD9Oig0gqm"
      },
      "source": [
        "* Change the default parameters of the auto_arima function to get faster results. Eploit the fact that you know the \"true\" data generating process."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "rOlZ7ufT0gqm"
      },
      "outputs": [],
      "source": [
        "model = auto_arima(\n",
        "    y,\n",
        "    d=?,\n",
        "    start_p=?,\n",
        "    max_p=?,\n",
        "    start_q=?,\n",
        "    max_q=?,\n",
        "    D=?,\n",
        "    start_P=?,\n",
        "    max_P=?,\n",
        "    start_Q=?,\n",
        "    max_Q=?,\n",
        "    seasonal=?,\n",
        "    information_criterion=\"bic\",\n",
        ")\n",
        "\n",
        "result_table(model)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "b-sJZqhqTctH"
      },
      "outputs": [],
      "source": [
        "model.summary()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ePVK6Zdh0gqm"
      },
      "source": [
        "## Real-world time series"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "nIt5lBaK0gqm"
      },
      "outputs": [],
      "source": [
        "# Read the data file\n",
        "DF = pd.read_csv(\"../data/public/train.csv\", parse_dates=True)\n",
        "DF[\"Day\"] = pd.to_datetime(DF[\"Day\"], format=\"%Y-%m-%d\")\n",
        "DF.set_index(\"Day\", inplace=True)\n",
        "DF = DF.asfreq(\"D\")\n",
        "DF.fillna(method=\"backfill\", inplace=True)\n",
        "\n",
        "DF"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "1KAUQ5nQ0gqn"
      },
      "outputs": [],
      "source": [
        "# Select the series to consider\n",
        "id_series_all = [\"s001\", \"s002\", \"s003\"]\n",
        "DF_all = DF[id_series_all].copy()\n",
        "\n",
        "DF_all"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "iL3xzLCT0gqn"
      },
      "outputs": [],
      "source": [
        "HORIZON = 7 * 2\n",
        "DF_train = DF_all[:-HORIZON]\n",
        "DF_test = DF_all[-HORIZON:]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xoVtupmt0gqn"
      },
      "source": [
        "* Use auto.arima to compute forecasts for all series"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "bzWs-f7J0gqn"
      },
      "outputs": [],
      "source": [
        "from pmdarima.arima import auto_arima\n",
        "\n",
        "fcts_arima_list = list()\n",
        "\n",
        "for id_series in id_series_all:\n",
        "    print(\"======\", id_series, \"======\")\n",
        "    y = DF_train[id_series]\n",
        "    model = auto_arima(\n",
        "        y,\n",
        "        d=?,\n",
        "        start_p=?,\n",
        "        max_p=?,\n",
        "        start_q=?,\n",
        "        max_q=?,\n",
        "        D=?,\n",
        "        start_P=?,\n",
        "        max_P=?,\n",
        "        start_Q=?,\n",
        "        max_Q=?,\n",
        "        m=?,\n",
        "        trace=True,\n",
        "    )\n",
        "    f_arima = model.predict(HORIZON)\n",
        "    f_arima.name = id_series\n",
        "    fcts_arima_list.append(f_arima)\n",
        "\n",
        "fcts_arima = pd.concat(fcts_arima_list, axis=1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZNoJpAVz0gqn"
      },
      "source": [
        "* Compute naive forecasts for all series"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ixWVr6nP0gqo"
      },
      "outputs": [],
      "source": [
        "fcts_mean_list = list()\n",
        "fcts_naive_list = list()\n",
        "fcts_snaive_list = list()\n",
        "\n",
        "for series in id_series_all:\n",
        "    series_train = DF_train[series]\n",
        "    series_test = DF_test[series]\n",
        "\n",
        "    T = len(series_train)\n",
        "\n",
        "    ## Mean\n",
        "    meanf = series_train.mean()\n",
        "    f_mean = pd.Series([meanf for h in range(0, HORIZON)], index=series_test.index)\n",
        "    f_mean.name = series\n",
        "\n",
        "    fcts_mean_list.append(f_mean)\n",
        "\n",
        "    ## Naive\n",
        "    f_naive = series_train[-1]\n",
        "    f_naive = pd.Series([f_naive for h in range(0, HORIZON)], index=series_test.index)\n",
        "    f_naive.name = series\n",
        "\n",
        "    fcts_naive_list.append(f_naive)\n",
        "\n",
        "    ## Seasonal naive\n",
        "    # period = 7\n",
        "    # f_snaive = [series_train[T + h - period * ((HORIZON -1)//period + 1)] for h in range(0, HORIZON) ]\n",
        "    f_snaive = [series_train[-HORIZON + h] for h in range(0, HORIZON)]\n",
        "    f_snaive = pd.Series(f_snaive, index=series_test.index)\n",
        "    f_snaive.name = series\n",
        "\n",
        "    fcts_snaive_list.append(f_snaive)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "n_KHI_CA0gqo"
      },
      "outputs": [],
      "source": [
        "fcts_mean = pd.concat(fcts_mean_list, axis=1)\n",
        "fcts_naive = pd.concat(fcts_naive_list, axis=1)\n",
        "fcts_snaive = pd.concat(fcts_snaive_list, axis=1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "86mLF1P-0gqo"
      },
      "source": [
        "* Compute the Symmetric Mean Absolute Percentage Error (SMAPE) for all series."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "9hKcRLax0gqo"
      },
      "outputs": [],
      "source": [
        "def smape(y_true, y_pred):\n",
        "    assert (y_pred >= 0).all().all()\n",
        "    denominator = (y_true + y_pred) / 200.0\n",
        "    SAPE = np.abs(y_true - y_pred) / denominator\n",
        "    SAPE[denominator == 0] = 0.0\n",
        "    return SAPE.mean().mean()\n",
        "\n",
        "\n",
        "print(smape(DF_test, fcts_mean))\n",
        "print(smape(DF_test, fcts_naive))\n",
        "print(smape(DF_test, fcts_snaive))\n",
        "print(smape(DF_test, fcts_arima))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "durZWcgU0gqo"
      },
      "outputs": [],
      "source": [
        "id_series = \"s001\"\n",
        "fcts_snaive[id_series].plot(label=\"snaive\")\n",
        "DF_test[id_series].plot(label=\"true\")\n",
        "fcts_arima[id_series].plot(label=\"arima\")\n",
        "plt.legend()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Gs50dttz0gqo"
      },
      "source": [
        "# Neural network forecast\n",
        "\n",
        "In the following, we will consider two neural network architectures for forecasting. Your task is to play with all the hyperparameters to obtain the best out-of-sample forecasts, i.e. on the test set.\n",
        "\n",
        "Some important hyperparameters include: n_simul (size of the dataset), LAG (the number of lagged values), LATENT_DIM (the number of units in the layer), BATCH_SIZE (number of samples per mini-batch), EPOCHS (the number of epochs), the optimizer and the early stop strategy."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "OOWAWxp10gqo"
      },
      "outputs": [],
      "source": [
        "from pmdarima.arima import auto_arima\n",
        "import keras\n",
        "from keras.callbacks import EarlyStopping, ModelCheckpoint\n",
        "\n",
        "from main.utils.utils_methods import embed_data, plot_learning_curves\n",
        "from main.utils.utils import mse, mae, mape, smape\n",
        "\n",
        "%load_ext autoreload\n",
        "%autoreload 2"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6McRK5sp0gqo"
      },
      "source": [
        "* We will simulate a time series from a nonlinear stochastic process:\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "IragzmJz0gqo"
      },
      "outputs": [],
      "source": [
        "n_simul = 1000\n",
        "n_burn = 100\n",
        "n = n_simul + n_burn\n",
        "noise = np.random.normal(size=n)\n",
        "\n",
        "y = np.zeros(n)\n",
        "y[0] = 0\n",
        "y[1] = 0\n",
        "for t in range(2, n):\n",
        "    y[t] = (\n",
        "        0.3 * y[t - 1]\n",
        "        + 0.6 * y[t - 2]\n",
        "        + (0.1 - 0.9 * y[t - 1] + 0.8 * y[t - 2]) * (1 / (1 + np.exp(-10 * y[t - 1])))\n",
        "        + noise[t]\n",
        "    )\n",
        "\n",
        "data = pd.DataFrame(y[n_burn:], columns=[\"series\"])\n",
        "plt.plot(data)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3pquchMO0gqp"
      },
      "source": [
        "Choose which loss function you want to experiment with. It is used later in the code to fit and evaluate a neural network model."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "dGXoB0ea0gqp"
      },
      "outputs": [],
      "source": [
        "# Loss function to be used to optimize the model parameters\n",
        "loss_fct = \"mse\"  # 'mae'\n",
        "\n",
        "# Accuracy measure to be used to evaluate test predictions.\n",
        "accuracy_measure = mse  # mae # mape # smape"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "noN3tVMX0gqp"
      },
      "outputs": [],
      "source": [
        "# The forecast horizon\n",
        "HORIZON = 3\n",
        "\n",
        "# The number of lagged values.\n",
        "LAG = 4\n",
        "\n",
        "# Data split\n",
        "n = len(data)\n",
        "n_train = int(0.6 * n)\n",
        "n_valid = int(0.2 * n)\n",
        "n_learn = n_train + n_valid\n",
        "\n",
        "train = data[:n_train]\n",
        "valid = data[n_train:n_learn]\n",
        "test = data[n_learn:]\n",
        "\n",
        "# From time series to input-output data (also called time series embedding)\n",
        "(\n",
        "    train_inputs,\n",
        "    valid_inputs,\n",
        "    test_inputs,\n",
        "    X_train,\n",
        "    y_train,\n",
        "    X_valid,\n",
        "    y_valid,\n",
        "    X_test,\n",
        "    y_test,\n",
        ") = embed_data(train, valid, test, HORIZON, LAG, freq=None, variable=\"series\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "23A89AkIEFt2"
      },
      "outputs": [],
      "source": [
        "display(X_train.head())\n",
        "display(y_train.head())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "F1Z5dWXN0gqp"
      },
      "source": [
        "# Multioutput MLP"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "X81PSGnDX3jR"
      },
      "outputs": [],
      "source": [
        "from keras.models import Sequential\n",
        "from keras.callbacks import EarlyStopping, ModelCheckpoint\n",
        "from keras.layers import Dense\n",
        "\n",
        "# Simple MLP with 1 hidden layer\n",
        "def mlp_multioutput(X_train, y_train, X_valid, y_valid, LATENT_DIM, BATCH_SIZE, EPOCHS, LAG, HORIZON, \n",
        "                    loss, optimizer, earlystop, best_val, verbose):\n",
        "  \n",
        "    model = Sequential()\n",
        "    model.add(Dense(LATENT_DIM, activation=\"relu\", input_shape=(LAG,)))\n",
        "    model.add(Dense(HORIZON))\n",
        "    model.compile(optimizer=optimizer, loss=loss)\n",
        "    \n",
        "    history = model.fit(\n",
        "        X_train,\n",
        "        y_train,\n",
        "        batch_size=BATCH_SIZE,\n",
        "        epochs=EPOCHS,\n",
        "        validation_data=(X_valid, y_valid),\n",
        "        callbacks=[earlystop, best_val],\n",
        "        verbose=verbose\n",
        "    )\n",
        "    return model, history"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "jBvIJFAn0gqp"
      },
      "outputs": [],
      "source": [
        "#########################\n",
        "file_header = \"model_\" + \"mlp_multioutput\"\n",
        "verbose = 0\n",
        "\n",
        "optimizer_adam = keras.optimizers.Adam(learning_rate=0.01)\n",
        "earlystop = EarlyStopping(monitor=\"val_loss\", min_delta=0, patience=50)\n",
        "\n",
        "LATENT_DIM = 5  # 50   # number of units in the RNN layer\n",
        "BATCH_SIZE = 32  # number of samples per mini-batch\n",
        "EPOCHS = 200  # maximum number of times the training algorithm will cycle through all samples\n",
        "loss = loss_fct\n",
        "\n",
        "best_val = ModelCheckpoint(\n",
        "    \"../work/\" + file_header + \"_{epoch:02d}.h5\",\n",
        "    save_best_only=True,\n",
        "    mode=\"min\",\n",
        "    save_freq=\"epoch\",\n",
        "    monitor=\"val_loss\",\n",
        ")\n",
        "#########################\n",
        "\n",
        "model_mlp_multioutput, history_mlp_multioutput = mlp_multioutput(\n",
        "    X_train,\n",
        "    y_train,\n",
        "    X_valid,\n",
        "    y_valid,\n",
        "    LATENT_DIM=LATENT_DIM,\n",
        "    BATCH_SIZE=BATCH_SIZE,\n",
        "    EPOCHS=EPOCHS,\n",
        "    LAG=LAG,\n",
        "    HORIZON=HORIZON,\n",
        "    loss=loss,\n",
        "    optimizer=optimizer_adam,\n",
        "    earlystop=earlystop,\n",
        "    best_val=best_val,\n",
        "    verbose=verbose,\n",
        ")\n",
        "\n",
        "plot_learning_curves(history_mlp_multioutput)\n",
        "\n",
        "best_epoch = np.argmin(np.array(history_mlp_multioutput.history[\"val_loss\"])) + 1\n",
        "print(\"Best epoch:\", best_epoch)\n",
        "filepath = \"../work/\" + file_header + \"_{:02d}.h5\"\n",
        "model_mlp_multioutput.load_weights(filepath.format(best_epoch))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "A0P-Sqy00gqq"
      },
      "source": [
        "## Recursive MLP\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "KK1avMBj0gqq"
      },
      "outputs": [],
      "source": [
        "#########################\n",
        "file_header = \"model_\" + \"mlp_recursive\"\n",
        "verbose = 0\n",
        "\n",
        "optimizer_adam = keras.optimizers.Adam(learning_rate=0.01)\n",
        "earlystop = EarlyStopping(monitor=\"val_loss\", min_delta=0, patience=50)\n",
        "\n",
        "LATENT_DIM = 5  # number of units in the RNN layer\n",
        "BATCH_SIZE = 32  # number of samples per mini-batch\n",
        "EPOCHS = 200  # maximum number of times the training algorithm will cycle through all samples\n",
        "loss = loss_fct\n",
        "\n",
        "best_val = ModelCheckpoint(\n",
        "    \"../work/\" + file_header + \"_{epoch:02d}.h5\",\n",
        "    save_best_only=True,\n",
        "    mode=\"min\",\n",
        "    save_freq=\"epoch\",\n",
        "    monitor=\"val_loss\",\n",
        ")\n",
        "#########################\n",
        "\n",
        "(\n",
        "    _,\n",
        "    _,\n",
        "    _,\n",
        "    X_train_onestep,\n",
        "    y_train_onestep,\n",
        "    X_valid_onestep,\n",
        "    y_valid_onestep,\n",
        "    _,\n",
        "    _,\n",
        ") = embed_data(train, valid, test, 1, LAG, freq=None, variable=\"series\")\n",
        "\n",
        "# The recursive MLP is just a multioutput MLP with 1 output.\n",
        "# However, the predictions given by the multioutput MLP and recursive MLP \n",
        "# are not created in the same way.\n",
        "model_mlp_recursive, history_mlp_recursive = mlp_multioutput(\n",
        "    X_train_onestep,\n",
        "    y_train_onestep,\n",
        "    X_valid_onestep,\n",
        "    y_valid_onestep,\n",
        "    LATENT_DIM=LATENT_DIM,\n",
        "    BATCH_SIZE=BATCH_SIZE,\n",
        "    EPOCHS=EPOCHS,\n",
        "    LAG=LAG,\n",
        "    HORIZON=1,\n",
        "    loss=loss,\n",
        "    optimizer=optimizer_adam,\n",
        "    earlystop=earlystop,\n",
        "    best_val=best_val,\n",
        "    verbose=verbose,\n",
        ")\n",
        "plot_learning_curves(history_mlp_recursive)\n",
        "\n",
        "best_epoch = np.argmin(np.array(history_mlp_recursive.history[\"val_loss\"])) + 1\n",
        "print(\"Best epoch:\", best_epoch)\n",
        "filepath = \"../work/\" + file_header + \"_{:02d}.h5\"\n",
        "model_mlp_recursive.load_weights(filepath.format(best_epoch))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "K-DFLG5U0gqq"
      },
      "source": [
        "# Naive forecasts"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "5TQiDULj0gqq"
      },
      "outputs": [],
      "source": [
        "# len(X_test.values[:, -1])\n",
        "predictions_naive = np.tile(X_test.values[:, -1], (HORIZON, 1)).T\n",
        "predictions_naive = pd.DataFrame(\n",
        "    predictions_naive, columns=[f\"t+{t}\" for t in range(1, HORIZON + 1)]\n",
        ")\n",
        "predictions_naive"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5k5qG7vf0gqq"
      },
      "source": [
        "# MLP forecasts"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "iFbfh-qseOdk"
      },
      "outputs": [],
      "source": [
        "predictions_mlp_multioutput = model_mlp_multioutput.predict(X_test)\n",
        "predictions_mlp_multioutput = pd.DataFrame(\n",
        "    predictions_mlp_multioutput, columns=[f\"t+{t}\" for t in range(1, HORIZON + 1)]\n",
        ")\n",
        "\n",
        "predictions_mlp_multioutput"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "vduiEE9xfQ65"
      },
      "outputs": [],
      "source": [
        "X_test"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "D2_HwzZ-0gqr"
      },
      "outputs": [],
      "source": [
        "for h in range(HORIZON):\n",
        "    pred = model_mlp_recursive.predict(X_test)\n",
        "    # `predictions_mlp_recursive` contains the predictions\n",
        "    if h == 0:\n",
        "        predictions_mlp_recursive = pred\n",
        "    else:\n",
        "        predictions_mlp_recursive = np.hstack((predictions_mlp_recursive, pred))\n",
        "    # `X_test` is updated at each step\n",
        "    X_test = pd.DataFrame(\n",
        "        np.hstack((X_test.to_numpy()[:, 1:], pred)),\n",
        "        index=X_test.index,\n",
        "        columns=X_test.columns,\n",
        "    )\n",
        "\n",
        "predictions_mlp_recursive = pd.DataFrame(\n",
        "    predictions_mlp_recursive, columns=[f\"t+{t}\" for t in range(1, HORIZON + 1)]\n",
        ")\n",
        "\n",
        "predictions_mlp_recursive"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "m9Ic1AWYeVA1"
      },
      "outputs": [],
      "source": [
        "predictions_combination = (predictions_mlp_multioutput + predictions_mlp_recursive) / 2"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3-8BKotM0gqr"
      },
      "source": [
        "# ARIMA forecasts"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "AVDWYLca0gqr"
      },
      "outputs": [],
      "source": [
        "model = auto_arima(data[:n_learn])\n",
        "fcts_list = []\n",
        "for i in np.arange(len(y_test)):\n",
        "    pred = model.fit_predict(data[LAG + i : n_learn + LAG + i], n_periods=HORIZON)\n",
        "    fcts_list.append(pred.to_numpy()[np.newaxis])\n",
        "\n",
        "predictions_arima = pd.DataFrame(\n",
        "    np.concatenate(fcts_list), columns=[\"t+\" + str(t) for t in range(1, HORIZON + 1)]\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aAdHydZR0gqr"
      },
      "source": [
        "# Forecast accuracy"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "A6tT7-pD0gqr"
      },
      "outputs": [],
      "source": [
        "true_values = pd.DataFrame(\n",
        "    test_inputs[\"target\"], columns=[\"t+\" + str(t) for t in range(1, HORIZON + 1)]\n",
        ")\n",
        "\n",
        "predictions = {\n",
        "    'naive': predictions_naive,\n",
        "    'mlp_multioutput': predictions_mlp_multioutput,\n",
        "    'mlp_recursive': predictions_mlp_recursive,\n",
        "    'combination': predictions_combination,\n",
        "    'arima': predictions_arima,\n",
        "}\n",
        "\n",
        "results = {}\n",
        "for model_name, prediction in predictions.items():\n",
        "    results[model_name] = []\n",
        "    for h in range(1, HORIZON + 1):\n",
        "        time_horizon = \"t+\" + str(h)\n",
        "        results[model_name].append(\n",
        "            accuracy_measure(true_values[time_horizon], prediction[time_horizon])\n",
        "        )\n",
        "\n",
        "pd.DataFrame(results).mean().to_frame()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "UfNEuXmCZA9a"
      },
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3.8.9 64-bit",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.8.9"
    },
    "orig_nbformat": 2,
    "vscode": {
      "interpreter": {
        "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}