In questo post indagheremo la possibilità di utilizzare tecniche di machine learning e in particolare algoritmi di classificazione per l'analisi dei mercati azionari e il raggruppamento di titoli in base alla loro esposizione alle componenti principali dei mercati, che abbiamo già visto in precedenza. Utilizziamo a tal fine i primi venti titoli che compongono l'indice FTSEMIB dei quali scarichiamo da Yahoo! Finance le serie storiche degli ultimi due anni, trasformandole infine in rendimenti giornalieri.

import pandas as pd
import numpy as np
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
import altair as alt
import yfinance as yf

tickers = [
    'ENEL.MI',
    'ISP.MI',
    'STLA.MI',
    'ENI.MI',
    'UCG.MI',
    'G.MI',
    'RACE.MI',
    'STM.MI',
    'CNHI.MI',
    'MONC.MI',
    'SRG.MI',
    'FBK.MI',
    'TRN.MI',
    'NEXI.MI',
    'PRY.MI',
    'ATL.MI',
    'MB.MI',
    'EXO.MI',
    'CPR.MI',
    'AMP.MI'
]

data = yf.download(tickers=tickers, period='2y')
prices = data.xs('Adj Close', axis=1, level=0)
returns = prices.pct_change()
returns.head()

[*********************100%***********************]  20 of 20 completed
AMP.MI ATL.MI CNHI.MI CPR.MI ENEL.MI ENI.MI EXO.MI FBK.MI G.MI ISP.MI MB.MI MONC.MI NEXI.MI PRY.MI RACE.MI SRG.MI STLA.MI STM.MI TRN.MI UCG.MI
Date
2019-06-28 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2019-07-01 0.001945 -0.032300 0.013972 0.012188 -0.002280 0.002328 0.007792 0.005097 0.004227 -0.007651 -0.012795 0.010638 0.003532 0.007714 0.008754 -0.004805 0.006685 0.041987 -0.005714 -0.014964
2019-07-02 0.000000 0.037889 0.000437 0.014335 0.021874 -0.009699 0.008376 0.005071 0.008118 0.009317 0.011397 -0.002105 -0.002419 0.000546 0.013884 0.026437 -0.003725 -0.009535 0.035201 -0.010690
2019-07-03 0.021359 0.024772 0.026891 0.020916 0.023642 0.004828 0.021406 0.046922 0.017895 0.050822 0.026072 0.015295 0.007496 -0.000546 0.017460 0.020829 0.028126 -0.005901 0.018390 0.055735
2019-07-04 -0.007605 0.005937 0.010219 -0.019380 0.006086 0.005354 0.007194 0.031807 0.002051 0.022819 0.019810 0.002857 -0.007222 0.004647 -0.002692 -0.006362 0.003004 0.008435 -0.004429 0.049560

Calcoliamo poi la matrice di correlazione tra i rendimenti dei titoli e le componenti principali di quest'ultima, utilizzando per l'analisi le prime sei componenti, che spiegano circa l'80% - 85% della varianza.

corr_mat = returns.corr()

n_components = 6
pca = PCA(n_components=n_components)
pc = pca.fit_transform(corr_mat)
pc_df = pd.DataFrame(
    pc,
    returns.columns,
    [f'PC{c}' for c in range(n_components)]
)
pc_df.index.name = 'Titolo'
pc_df.head()

PC0 PC1 PC2 PC3 PC4 PC5
Titolo
AMP.MI 0.690079 -0.278930 0.001941 -0.165962 -0.210654 -0.184466
ATL.MI 0.122440 -0.401362 0.664491 0.236238 0.019327 0.070195
CNHI.MI -0.373885 -0.201208 0.139796 -0.158293 0.091624 0.253066
CPR.MI 0.634638 0.050509 0.136282 -0.330863 -0.010690 -0.134976
ENEL.MI 0.247968 0.389905 -0.023243 0.140282 -0.121522 0.048196

var_df = pd.DataFrame(
    pca.explained_variance_ratio_,
    pc_df.columns,
    ['Varianza cumulata']
)
var_df.index.name = 'PC'
alt.Chart(
    var_df.cumsum().reset_index()
).mark_bar().encode(
    y='PC:N',
    x='Varianza cumulata:Q',
    tooltip=['PC', 'Varianza cumulata']
).properties(
    title='Varianza cumulata spiegata dalle componenti principali'
).interactive()

Possiamo analizzare quali sono le esposizioni dei singoli titoli alle componenti principali considerate - e in particolare alla prima, che spiega circa il 40% della varianza. Vediamo che le azioni si distribuiscono circa in due metà, con esposizioni negative e positive rispetto alla prima componente.

alt.Chart(
    pc_df.reset_index().melt(
        'Titolo',
        var_name='PC',
        value_name='Esposizione'
    )
).mark_bar().encode(
    y='Titolo:N',
    x='Esposizione:Q',
    color='PC',
    row='PC',
    tooltip=['Titolo', 'PC', 'Esposizione']
).properties(
    title='Esposizione dei titoli alle componenti principali'
).interactive()

Facendo il grafico dei valori delle componenti principali, possiamo osservare che la prima componente è molto più volatile delle altre - per evitare che questo fenomeno influenzi la classificazione dei titoli, standardizziamo) i rendimenti delle componenti principali.

pc_px = returns.dot(pc_df).add(1).cumprod()
pc_px.iloc[0, :] = 1
alt.Chart(
    pc_px.reset_index().melt(
        'Date',
        var_name='PC',
        value_name='Prezzo'
    )
).mark_line().encode(
    x='Date:T',
    y='Prezzo:Q',
    color='PC:N',
    tooltip=['Date', 'PC', 'Prezzo']
).properties(
    title='Andamento delle componenti principali'
).interactive()

def zscore(x):
    return (x - x.mean()) / x.std()

zpc_df = pc_df.apply(zscore)

Applichiamo a questo punto l'algoritmo K-means per classificare i titoli in K gruppi simili. Dovendo decidere ex-ante il numero K di gruppi in cui partizionare lo spazio dei titoli, utilizziamo due metodi di ricerca dell'ottimo applicandoli a valori di K compresi tra 2 e 10. Il primo metodo (inertia) fissa il punto di ottimo al "gomito" del grafico ottenuto - la linea più in alto nella figura sottostante. Il secondo metodo (silhouette) trova il numero di cluster ottimo in corrispondenza del massimo (locale) della funzione. In entrambi i casi, 7 sembra essere il valore ottimale per noi.

metrics_df = pd.DataFrame(
    index=range(2, 10),
    columns=['Inertia', 'Silhouette']
)
metrics_df.index.name = 'N_clusters'
for i in metrics_df.index:
    km = KMeans(n_clusters=i)
    km.fit_predict(zpc_df)
    metrics_df.loc[i, 'Inertia'] = km.inertia_
    metrics_df.loc[i, 'Silhouette'] = 100 * silhouette_score(
        zpc_df, 
        km.labels_, 
        metric='euclidean'
    )
metrics_df

Inertia Silhouette
N_clusters
2 96.27442 14.174693
3 79.545295 17.737655
4 64.363608 23.572529
5 48.211426 23.755179
6 33.114174 32.99203
7 22.553873 35.571988
8 18.496865 33.401054
9 15.065332 30.904044

alt.Chart(
    metrics_df.reset_index().melt(
        'N_clusters',
        var_name='Metric',
        value_name='Score'
    )
).mark_line().encode(
    x='N_clusters:O',
    y='Score:Q',
    color='Metric:N',
    row='Metric:N',
    tooltip=['N_clusters', 'Score']
).properties(
    title='Metriche per la scelta del numero ottimale di clusters'
).interactive()

Raggruppiamo quindi i nostri venti titoli in 7 cluster utilizzando l'algoritmo K-means. Essendo cluster multidimensionali, è complesso rappresentarli graficamente; tuttavia, è possibile costruire uno scatterplot che visualizzi i titoli raggruppati nei cluster sulla base delle prime due componenti principali. Possiamo notare nel grafico due cluster al margine sinistro che comprendono rispettivamente le grandi banche più Generali e Eni e i titoli collegati al gruppo Stellantis (ma anche Ferrari e STM, collocati in posizione più centrale). Altri due cluster facilmente identificabili sono nella parte alta del grafico (le aziende di pubblica utilità) e al margine destro (Campari e Amplifon, più Moncler, che si posiziona invece al centro del grafico). Infine gli ultimi tre cluster sono meno facilmente identificabili in base alle prime due componenti principali e comprendono rispettivamente Finecobank e Prysmian, la sola Atlantia e la sola Nexi.

n_clusters = 7
km_opt = KMeans(n_clusters=n_clusters)
clusters = km_opt.fit_predict(zpc_df)
cl_df = pc_df.copy()
cl_df['Cluster'] = clusters
cl_df.head()

PC PC0 PC1 PC2 PC3 PC4 PC5 Cluster
Titolo
AMP.MI 0.690079 -0.278930 0.001941 -0.165962 -0.210654 -0.184466 0
ATL.MI 0.122440 -0.401362 0.664491 0.236238 0.019327 0.070195 2
CNHI.MI -0.373885 -0.201208 0.139796 -0.158293 0.091624 0.253066 5
CPR.MI 0.634638 0.050509 0.136282 -0.330863 -0.010690 -0.134976 0
ENEL.MI 0.247968 0.389905 -0.023243 0.140282 -0.121522 0.048196 4

alt.Chart(
    cl_df.reset_index(),
).mark_circle().encode(
    x='PC0:Q',
    y='PC1:Q',
    color='Cluster:N',
    tooltip=['Titolo', 'Cluster']
).properties(
    title='Cluster rispetto alle prime due componenti principali'
).interactive()

Abbiamo così identificato gruppi di titoli omogenei in base a caratteristiche intrinseche del mercato, piuttosto che in base a elementi anagrafici o descrittivi come possono essere la nazione o il settore di appartenenza. Potremmo ad esempio costruire un portafoglio in cui ciascun cluster ottiene lo stesso peso e confrontarne l'andamento con un portafoglio equipesato in base ai titoli.

wgt_df = pd.DataFrame(0.20 / cl_df['Cluster'].value_counts()).rename(columns={'Cluster': 'Peso'})

pesi = pd.DataFrame(index=cl_df.index, columns=['Equal', 'Clustered'])
pesi.Equal = 0.05
pesi.Clustered = cl_df.Cluster.map(wgt_df.Peso)

test_df = pd.concat([
    returns.dot(pesi.Equal).add(1).cumprod(),
    returns.dot(pesi.Clustered).add(1).cumprod()
], axis=1).rename(columns={0: 'Equal weights', 1: 'Clustered weights'})
alt.Chart(
    test_df.reset_index().melt(
        'Date',
        var_name='Portafoglio',
        value_name='Prezzo'
    )
).mark_line().encode(
    x='Date:T',
    y='Prezzo:Q',
    color='Portafoglio:N',
    tooltip=['Date', 'Prezzo']
).properties(
    title='Backtest di due portafogli basati sui clusters'
).interactive()

Risulta che il portafoglio equipesato rispetto ai cluster ha ottenuto un risultato migliore negli ultimi due anni rispetto al portafoglio equipesato sui titoli, ma al prezzo di una maggiore volatilità. Naturalmente, non si tratta di una tecnica di ottimizzazione robusta, ma di un algoritmo che permette di individuare nuove dimensioni di analisi dei titoli e dei mercati che possono costituire la base per indagare più a fondo le motivazioni della similarità o dissimilarità di date azioni.