Prédiction de la volatilité sur les marchés US - CFM Challenge 2018

by Romain Ait Abdelmalek-Lomenech et Robert Dazi

Présentation du jeu de données

Introduction

Dans le cadre de la gestion de capital, le risque associé a un port-folio d'actif est une quantité primordiale dans l'équilibrage de ce dernier. Dans le cadre de stratégies plus aggressives, la volatilité de fin séance est une variable d'une grande importance pour des raisons de gestion du risque et d'optimisation des résultats. Capital Fund Management (CFM) a ainsi en 2018 proposé comme défi d'apprentissage automatique l'estimation de volatilité de cloture sur les marchés actions US. Il existe de nombreuses manières de calculer la volatilité d'un actif et CFM a imposé sa propre formule non divulguée, la volatilité dépendant de la modélisation que l'on considère la mieux adaptée pour un actif et une stratégie donnés.

Le jeu de données est composé, pour un ensemble d'actions (318 précisement) issues du marche américain de la volatilité et du signe du rendement (au sens de CFM) toutes les 5 minutes de 9h30 à 14h pour un jour donnée. A partir de ces données historiques de volatilités devait être prédit la volatilité moyenne des deux dernières heures avant clôture de la journée. Ainsi on dispose de séries temporelles pour un ensemble de date et de produits fixés (ces classes sont approximativement équi-distribuées), chacune de ces séries étant associée à une grandeur réelle positive à prédire.

Nous disposons, de plus, d'un second jeu de données test, où les cibles sont absentes, et dont la précision des prédictions est obtenue en soumettant nos résultats à CFM.

De manière surprenante, le choix de la métrique est la MAPE, qui est connue pour son instabilité et qui est malheureusement peu supportée par les bibliothèques standards de Machine Learning.

De l'aveu de l'organisateur lui-même, faire sensiblement mieux qu'une régression linéaire aurait très difficile. Ce constat est partagé à la suite de notre étude. Étant un challenge déjà terminé, le meilleures score est public et s´élève à 20.7%, obtenu semble-t-il par stacking de réseaux de neuronnes.

Ce challenge nous a permis d'améliorer notre méthodologie, mesurer la difficulté d'optimiser un modèle qui marche déjà bien, nous a réservé des surprises en termes de compromis biais-variance (rajouter des colonnes peut sensiblement détériorer un modèle) et nous a permis de s'essayer à des bibliothèques de Machine Learning montantes.

En raison de la taille de nos données, et comme nous le verrons plus loin du nombre de modèles à générer et optimiser, nous avons choisi de réaliser notre projet en python. Nous utilisons ainsi la bibliothèque H2o, qui couplée au cluster de calcul de l'université Paris-Sud, permet d'optimiser nos calculs et de les paralléliser

import io

import numpy as np
import pandas as pd
import dask.dataframe as dd
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns

from ipypublish import nb_setup

pd = nb_setup.setup_pandas(escape_latex=False)

from sklearn.linear_model import LinearRegression, Lars
from mlinsights.mlmodel import QuantileLinearRegression
from sklearn.preprocessing import OneHotEncoder

sns.set()
from submit import submit, submit_dataframe
# Loading the raw training data

train_raw = dd.read_csv("preprocess/training_without_imputation/*.part").compute(scheduler="processes")

test_raw = pd.read_csv("data/testing_input.csv", sep=";")
vol_cols = train_raw.columns[train_raw.columns.str.contains("volatility")]

Description succinte des données

Les données contiennent trois groupes de variables prédictrices:

  • volatilités: 54 colonnes qui represente les volatilités prises toutes les 5 minutes de 9h30 inclus à 14h exclus.
  • signe du rendement: 54 colonnes avec le même échantillonage que précedemment.
  • contexte: 3 colonnes identifiant du produit, identifiant de la date, identifiant unique de la série temporelle.

Les identifiants associés aux dates ont été mis au hasard sans considérer leur caractère séquentiel. Les identifiants des produits sont anonymisés.

Le nombre de date de quotation s'élève à 2117 jours, sachant qu'il y a 253 jours de quotation par an, on obtient un peu plus de 8 années de quotation.

Le fait que l'on ne puisse pas considérer les dates de manière séquentielle est handicapant. On ne peut ansi rajouter en terme de colonnes les volatilités du jour qui précéde ou une moyenne des volatilités sur la semaine. Il s'agit cependant d'un choix logique, étant donné qu'il serait possible d'utiliser la volatilité d'ouverture de séance pour aider à la prédiction de la volatilité de fin de séance de la veille, ce qui réduirait considérablement l'intérêt réel des prédictions.

Visualisation

Exemple d'une série temporelle obtenue

Nous en affichons ici trois et nous remarquons des comportements très différents en termes de pics de volatilités. Ainsi la première série a un pic en début de séance, tandis cela n'apparait qu'en fin d'après-midi pour la seconde, tandis que les pics de la troisième ont lieu en début et milieu de matinée.

De manière générale, les observations empiriques des marchés financiers décrites dans plusieurs articles permettent de constater une forme de "U" au cours de la journée, avec une volatilité forte en début de session, plus faible en milieu de journée, avec une remonté en fin de session, sans toutefois généralement atteindre la volatilité de début de session.

for ts_row in [57, 1999, 500]:
    ts_example = pd.DataFrame({
        "volatility": train_raw.iloc[ts_row, 3:57].values,
        "return": train_raw.iloc[ts_row, 57:111].values,
        "time": train_raw.columns[3:57].str[-8:],
    })

    plt.figure(figsize=(30, 10))
    sns.barplot("time", "volatility", hue="return", palette="RdBu", dodge=False, data=ts_example)
    plt.suptitle(f"Serie ID: {train_raw['ID'].iloc[ts_row]}, TARGET: {train_raw['TARGET'].iloc[ts_row]:.4}")
    plt.xticks(rotation=90);

En étudiant la distribution des volatilités conditionnellement à l'heure de la journée, on remarque ainsi, comme espéré, une très forte volatilité moyenne en début de séance.

L'affichage des outliers a été désactivé dans ce boxplot par soucis de visibilité. Les points triangulaires verts repésentent la moyenne de volatilité pour l'horaire considéré.

"""
ts_mean = train_raw.mean()
ts_mean_frame = pd.DataFrame({
    "volatility": ts_mean.iloc[3:57].values,
    "return": ts_mean.iloc[57:111].values,
    "time": ts_mean.index[3:57].str[-8:],
    })

plt.figure(figsize=(30, 10))
sns.barplot("time", "volatility", data=ts_mean_frame)
plt.suptitle(f"Serie ID: {train_raw['ID'].iloc[ts_row]}, TARGET: {train_raw['TARGET'].iloc[ts_row]:.4}")
plt.xticks(rotation=90);
"""
volatilities_per_time_span = train_raw[train_raw.columns[train_raw.columns.str.contains("volatility")]].melt()
volatilities_per_time_span["time"] = volatilities_per_time_span["variable"].str[10:].astype("category")

volatilities_per_time_span = volatilities_per_time_span.sort_values("time")

plt.figure(figsize=(30, 10))
sns.boxplot(x="time",
              y="value",
              order=volatilities_per_time_span["time"].drop_duplicates(),
              showmeans=True,
              showfliers = False,
              color="grey",
              data=volatilities_per_time_span)
plt.xticks(rotation=90);

Distribution des volatilités

Nous constatons visuellement que la log-volatilité sur les deux dernières heures avant fermeture suit presque une lois de student. Cela est en adéquation avec les observations empiriques du marché financier, où la log-volatilité bien que proche de la loi gaussienne affiche cependant un phénomène de fat-tail.

plt.figure(figsize=(20, 10))
sns.distplot(np.log(train_raw["TARGET"]));

De même, si on regarde les log-volatilités du matin, on obtient des distributions identiques en enlevant les sur-représentations en 0.

plt.figure(figsize=(20, 10))
sns.distplot(np.log(train_raw[vol_cols].melt().query("value > 0")["value"].dropna()));

Nous pouvons constater que la volatilité est nulle pour 4% des features :

train_raw[vol_cols].melt()["value"].eq(0).agg(["sum", "mean"]).to_frame()
value
sum 1.565805e+06
mean 4.556938e-02

Notons aussi la proportion de valeurs manquantes :

train_raw[vol_cols].melt()["value"].isna().agg(["sum", "mean"]).to_frame()
value
sum 419817.000000
mean 0.012218

Malheureusement, il est particulièrement difficile d'utiliser cet observation afin d'améliorer un modèle minimisant la MAPE.

Données Manquantes

Un des problèmes majeurs relatif aux données est la présence de nombreuses valeurs manquantes, aussi bien au niveau des features de chaque actif, que de la répartion des produits par date. L'absence de features s'explique pour deux raisons principales :

  • Instabilité du cours qui oblige la suspension de la quotation, mais cela reste particulièrement rare.
  • Aucun mouvement du cours, ce qui fait qu'il n'y a pas de mise a jour du prix du titre, ce qui indiquerait la présence de produits illiquides.

Pour l'imputation de ces données, nous avons adopté une interpolation linéaire pour les volatilités. Concernant les retours manquants, considérant la difficulté de les prédire de part l'aspect tout à fait aléatoire des retours des prix de marchés décrit dans la littérature, nous choisissons de fixer leur valeur à 0.

Un autre phénomène remarquable est présent au sein des données: pour certaines dates, certains produits sont manquants.

La encore, nous pouvons le justifier par:

  • action non encore quotée à la date donnée, ou retirée des marchés électroniques.
  • produit totalement illiquide à une date donnée.
  • volatilité nulle sur les deux heures qui impliquerait une inconsistence avec la métrique

Lien entre date et produit manquants

Afin de confirmer l'hypothese d'une très forte corrélation entre un ensemble de date et un ensemble de produits, nous calculons une matrice de presence entre produit et date et clusterisons cette matrice.

On remarque ci-bas en effet que sur l'ensemble d'entrainement les produits manquants sont concentrés sur un ensemble de date particulières.

cross_product_date_train = pd.crosstab(train_raw["product_id"], train_raw["date"])
sns.clustermap(cross_product_date_train,
              figsize=(20, 20),
              xticklabels=True,
             yticklabels=True,
              cmap="viridis");

Volatilités manquantes

Le graphe suivant a pour vocation de confirmer qu'en majorité les produits avec des proportions de valeurs manquantes élevés pour un jour donné sont ceux qui sont très peu liquide ce même jour.

sns.jointplot(train_raw[vol_cols].median(axis=1),
              train_raw[vol_cols].isna().mean(axis=1),
              xlim=(-0.1, 0.75),
              alpha=0.1,
              height=20);

Nous clusterisons la matrice qui à chaque produit et date associe la proportion de predicteurs manquants.

missing_prop_product_cross_date = pd.concat([train_raw[vol_cols].isna().mean(axis=1).rename("missing_prop"), train_raw[["product_id", "date"]]], axis=1).pivot_table(index="product_id", columns="date", values="missing_prop")

sns.clustermap(missing_prop_product_cross_date.fillna(0),
               figsize=(20, 20),
               xticklabels=True,
               yticklabels=True,
               cmap="viridis");