Impara a manipolare serie temporali con pandas e a condurre test di significatività tramite simulazione in Python, per analizzare la volatilità del mercato azionario.
Ottobre è storicamente il mese più volatile per le azioni, ma si tratta di un segnale persistente o solo rumore nei dati?
"Negli ultimi 32 anni, ottobre è stato in media il mese più volatile per l’S&P500 e dicembre il meno volatile".
In questo tutorial useremo Python per eseguire un’analisi completa e un test di questo fenomeno, per capire se è statisticamente significativo oppure no.
Useremo le funzionalità per serie temporali di pandas, come Resample, per trasformare i prezzi azionari grezzi in un formato che dia vita ai dati.
Scopriremo poi come condurre test d’ipotesi usando la simulazione in Python invece di formule macchinose.
Infine parleremo di un problema importante nell’analisi statistica, il bias da confronti multipli, vedremo come gestirlo e mostreremo visivamente con matplotlib l’effetto del “p-Hacking”.
Il nostro obiettivo:
Mostrare come usare Pandas per analizzare serie temporali
Capire come costruire un test d’ipotesi
Usare la simulazione con Python per eseguire test d’ipotesi
Mostrare l’importanza di considerare il bias da confronti multipli
I nostri dati:
Useremo dati giornalieri dell’S&P500 per questa analisi; in particolare, useremo i prezzi di chiusura giornalieri grezzi dal 1986 al 2018 (sorprendentemente difficili da trovare, quindi li ho resi pubblicamente disponibili).
L’ispirazione per questo post viene da Winton, che riprodurremo qui, seppur con 32 anni di dati contro i loro 87.
Lavorare i dati con Pandas
Per rispondere alla domanda se l’elevata volatilità osservata in alcuni mesi sia davvero significativa e quindi probabile che continui, dobbiamo trasformare i nostri 32 anni di dati sui prezzi in un formato che evidenzi il fenomeno che stiamo indagando.
Il formato scelto sarà il ranking medio mensile della volatilità (AMVR).
Il codice seguente mostra come portiamo i dati grezzi dei prezzi in questo formato. Iniziamo!
Per prima cosa le importazioni standard. (matplotlib.patches ci permette di controllare lo stile delle singole barre in un istogramma)
#standard importsimport numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
%matplotlib inline
Offerto da
Un trucco utile per far sì che i grafici occupino tutta la larghezza nei notebook Jupyter:
#resize charts to fit screen if using Jupyter Notebook
plt.rcParams['figure.figsize']=[15,5]Offerto da
Questa è opzionale: provala e guarda cosa succede!
Importiamo i dati usando il metodo read_csv, che prende un percorso, nel nostro caso un URL ai dati. Indichiamo anche di usare "date" come indice e di analizzare automaticamente le date dal testo fornito, al meglio possibile.
#Daily S&P500 data from 1986==>
url ="https://raw.githubusercontent.com/Patrick-David/Stocks_Significance_PHacking/master/spx.csv"
df = pd.read_csv(url,index_col='date', parse_dates=True)#view raw S&P500 data
df.head()Offerto da
close
date
1986-01-02
209.59
1986-01-03
210.88
1986-01-06
210.65
1986-01-07
213.80
1986-01-08
207.97
Questo ci dà i prezzi di chiusura grezzi e non rettificati dell’S&P500 (SPX). Ora dobbiamo convertire i prezzi grezzi in rendimenti giornalieri in percentuale. Per farlo, abbiamo due opzioni:
Potremmo prendere il logaritmo naturale dei prezzi. Questo dà un’approssimazione dei rendimenti giornalieri “veri”.
Possiamo usare il metodo pandas 'pct_change()' per calcolare direttamente la variazione percentuale giornaliera.
Per il nostro scopo useremo il metodo 2, dato che pandas può gestire istantaneamente il calcolo su un dataset di questa dimensione (oltre 8000 valori). Puliremo anche i dati eliminando il primo valore che diventa 'NaN' perché non c’è variazione di prezzo rispetto al giorno precedente. pct_change() accetta un parametro opzionale 'periods' che modifica lo shift temporale. Lo lasciamo all’impostazione predefinita pari a uno.
#To model returns we will use daily % change
daily_ret = df['close'].pct_change()#drop the 1st value - nan
daily_ret.dropna(inplace=True)#daily %change
daily_ret.head()Offerto da
date
1986-01-030.0061551986-01-06-0.0010911986-01-070.0149541986-01-08-0.0272681986-01-09-0.008944Name: close,dtype: float64
Offerto da
Il passo successivo è prendere queste variazioni giornaliere in percentuale e trasformarle in “volatilità mensile annualizzata”. La prima riga di codice sotto mostra che possiamo effettuare questa trasformazione in una sola riga. Questa è la potenza di Pandas! Ma scomponiamo ogni passaggio in sequenza.
1. Per ottenere mnthly_annu, usiamo innanzitutto il metodo 'resample' sui nostri rendimenti giornalieri. Resample ci permette di cambiare la frequenza dei periodi. Accetta una 'frequency offset string' o, in parole semplici, una lettera che corrisponde alla nuova frequenza desiderata. Includono:
B frequenza giorni lavorativi
C frequenza giorni lavorativi personalizzata
D frequenza giorni di calendario
W frequenza settimanale
M fine mese
Poiché vogliamo passare da rendimenti giornalieri a mensili, passiamo 'M' come parametro.
2. La seconda cosa è decidere come arrivare al nuovo valore mensile; ad esempio, potremmo sommare i valori giornalieri, moltiplicarli, ecc. Per la nostra analisi vogliamo una misura di volatilità e la deviazione standard funziona bene, quindi aggiungeremo std() ai nostri dati ricampionati per ottenere la volatilità mensile.
3. L’ultimo passo è annualizzare questo valore. Lo facciamo semplicemente moltiplicando per la radice quadrata di 12, con 12 che è il numero di periodi (mesi) in un anno. Questo ci dà i valori di volatilità mensile annualizzata che ci servono.
La buona pratica nell’analisi dei dati è visualizzarli man mano che si procede. Diamo quindi un’occhiata ai nostri dati di volatilità annualizzata.
Il grafico sotto mostra chiaramente i principali eventi di mercato, come il Black Monday e la Crisi Finanziaria del 2008. Il metodo matplotlib 'axvspan' ci consente di aggiungere linee verticali (axhspan è il corrispondente per le linee orizzontali); accetta 'xmin' e 'xmax' come parametri, che specificano l’ampiezza del rettangolo creato. Poiché il nostro asse x è l’indice datetime, possiamo semplicemente passare gli anni che vogliamo evidenziare. Il parametro alpha regola la trasparenza del blocco di colore per poter vedere il grafico sottostante.
L’attributo mpatches ci permette di creare una legenda personalizzata. Lo facciamo definendo prima la variabile 'labs' con colore, alpha e il testo desiderato, quindi passandola al parametro 'handles' in plt.legend per renderizzare la legenda.
#use pandas to resample returns per month and take Standard Dev as measure of Volatility#then annualize by multiplying by sqrt of number of periods (12)
mnthly_annu = daily_ret.resample('M').std()* np.sqrt(12)print(mnthly_annu.head())#we can see major market events show up in the volatility
plt.plot(mnthly_annu)
plt.axvspan('1987','1989',color='r',alpha=.5)
plt.axvspan('2008','2010',color='r',alpha=.5)
plt.title('Monthly Annualized vol - Black Monday and 2008 Financial Crisis highlighted')
labs = mpatches.Patch(color='red',alpha=.5, label="Black Monday & '08 Crash")
plt.legend(handles=[labs])Offerto da
date
1986-01-310.0333171986-02-280.0235851986-03-310.0279611986-04-300.0374261986-05-310.027412Freq:M,Name: close,dtype: float64
<matplotlib.legend.Legend at 0x280b1ee6908>Offerto da
Abbiamo quindi visto uno dei metodi potenti di pandas, resample. Ora usiamone un altro, 'groupby'. Dobbiamo passare dai valori di volatilità mensile annualizzata al nostro AMVR desiderato. Anche qui possiamo riuscirci in poche righe di codice.
1. Per prima cosa applichiamo 'groupby' ai valori di mnthly_annu. Groupby richiede un parametro che specifica come raggruppare i dati. Può essere una funzione o, come nel nostro caso, una Series. Passiamo 'mnthly_annu.index.year', che è semplicemente l’attributo anno dell’indice datetime di mnthly_annu. Questo raggrupperà i valori di volatilità mensile per ognuno dei 32 anni del nostro set di dati.
2. Successivamente, applichiamo il metodo 'rank', che ordina i dati in ordine crescente.
3. Infine, ripetiamo lo stesso e facciamo la media su tutti gli anni per ciascun mese. Questo ci dà i nostri valori AMVR finali!
#for each year rank each month based on volatility lowest=1 Highest=12
ranked = mnthly_annu.groupby(mnthly_annu.index.year).rank()#average the ranks over all years for each month
final = ranked.groupby(ranked.index.month).mean()
final.describe()Offerto da
count 12.000000
mean 6.450521
std 0.627458
min 5.21875025%6.03101350%6.49100475%6.704545
max 7.531250Name: close,dtype: float64
Offerto da
Questo ci dà il nostro ranking medio mensile della volatilità. Numericamente vediamo che il mese 10 (ottobre) è il più alto e il 12 (dicembre) è il più basso.
#the final average results over 32 years
final
Offerto da
date
16.81818226.66666736.57575847.30303056.60606166.03030376.03125085.87500096.406250107.531250116.343750125.218750Name: close,dtype: float64
Offerto da
Scegliere la giusta visualizzazione è importante perché vogliamo raccontare la storia dei dati con i grafici più efficaci. Per la nostra analisi vogliamo mostrare chiaramente la volatilità più alta e quella più bassa. Useremo quindi il grafico a barre di matplotlib, annotando e colorando di conseguenza per dare il massimo impatto.
1. Indicizzando b_plot (b_plot[9]) con la barra desiderata, possiamo impostarne il colore per evidenziare i valori massimo e minimo.
2. Per aggiungere i valori AMVR alla rispettiva barra, iteriamo su ciascun valore AMVR (variabile 'final') e arrotondiamo a 2 decimali. 'i' e 'v' rappresentano 'indice' e 'valore' per la variabile 'final', quindi 'i' andrà da 1 a 12 e 'value' dà il valore AMVR arrotondato a 2 decimali. Usiamo queste variabili in plt.text(). Il primo parametro di plt.text() è la posizione sull’asse x per l’etichetta, quindi passiamo 'i' per ogni barra e lo spostiamo di 0.8 per centrarla. Il secondo parametro è la stringa del valore di 'final'.
3. Per mostrare il valore medio, usiamo axhline e il parametro 'label' per aggiungere una legenda nell’angolo in alto a destra. Passiamo semplicemente il valore sull’asse y, cioè la media, e poi aggiungiamo lo stile.
Questo offre la visuale migliore del fenomeno che stiamo indagando. Si vede chiaramente che ottobre è stato il mese più volatile e dicembre il meno. Importante, si nota anche che dicembre è il valore più “estremo” in termini assoluti. Questo sarà importante per l’analisi nella sezione successiva.
#plot results for ranked s&p 500 volatility#clearly October has the highest AMVR#and December has the lowest#mean of 6.45 is plotted
b_plot = plt.bar(x=final.index,height=final)
b_plot[9].set_color('g')
b_plot[11].set_color('r')for i,v inenumerate(round(final,2)):
plt.text(i+.8,1,str(v), color='black', fontweight='bold')
plt.axhline(final.mean(),ls='--',color='k',label=round(final.mean(),2))
plt.title('Average Monthly Volatility Ranking S&P500 since 1986')
plt.legend()
plt.show()Offerto da
Questi sono i nostri dati; passiamo al test d’ipotesi…
Test d’ipotesi: qual è la domanda?
Il test d’ipotesi è una delle tecniche fondamentali della data science, eppure è tra le più intimidatorie e fraintese. La base di questa paura sta in come viene insegnato in Statistica 101, dove ci viene detto di:
"eseguire un t-test; è monodirezionale o bidirezionale? Scegli un test statistico adatto come il t-test di Welch, calcola i gradi di libertà, calcola il punteggio t, cerca il valore critico in tabella, confronta il valore critico con la statistica t…"
Comprensibilmente, questo porta a confusione su quale test eseguire e come eseguirlo. Tuttavia, tutte queste tecniche statistiche classiche per eseguire un test d’ipotesi sono state sviluppate quando avevamo pochissima potenza di calcolo ed erano semplicemente soluzioni analitiche chiuse per calcolare un p-value, tutto qui! Ma con l’ulteriore complicazione di dover scegliere la formula giusta per la situazione data, a causa delle loro ipotesi restrittive e talvolta opache.
Ma rallegriamoci!
C’è un modo migliore. La simulazione.
Per capire come la simulazione possa aiutarci, ricordiamo cos’è un test d’ipotesi:
Vogliamo testare “se l’effetto osservato nei nostri dati sia reale o se possa avvenire semplicemente per caso” e per eseguire questo test facciamo quanto segue:
Scegli un “test statistico” appropriato: è semplicemente un numero che misura l’effetto osservato. Nel nostro caso, sceglieremo la deviazione assoluta dell’AMVR dalla media.
Costruisci un’ipotesi nulla: è una versione dei dati in cui l’effetto osservato non è presente. Nel nostro caso, rimescoleremo ripetutamente le etichette dei dati (permutazione). La giustificazione è dettagliata sotto.
Calcola un p-value: è la probabilità di vedere l’effetto osservato tra i dati della nulla, in altre parole, per caso. Lo facciamo tramite simulazioni ripetute dei dati della nulla. Nel nostro caso, rimescoliamo più volte le etichette ‘date’ dei dati e contiamo semplicemente le occorrenze del nostro test statistico nelle diverse simulazioni.
Ecco il test d’ipotesi in 3 passi! Indipendentemente dal fenomeno che testiamo, la domanda è sempre la stessa: “l’effetto osservato è reale o dovuto al caso”
Il vero potere della simulazione è che dobbiamo rendere esplicite, attraverso il codice, le ipotesi del nostro modello. Mentre le tecniche classiche possono essere una “scatola nera” quanto alle loro ipotesi.
Esempio sotto: Il grafico a sinistra mostra i dati reali e l’effetto osservato con una certa probabilità (verde). Quello a destra è il nostro dato nullo simulato con la registrazione di quando l’effetto osservato è apparso per caso (rosso). Questa è la base del test d’ipotesi: qual è la probabilità di vedere l’effetto osservato sui dati della nulla.
La parte più critica del test d’ipotesi è chiarire che domanda stiamo cercando di rispondere. Nel nostro caso chiediamo:
“Il valore più estremo potrebbe verificarsi per caso?”
Il valore più estremo lo definiamo come la massima deviazione assoluta dell’AMVR dalla media. Questa domanda costituisce la nostra ipotesi nulla.
Nei nostri dati, il valore più estremo è quello di dicembre (1,23) non quello di ottobre (1,08), perché stiamo guardando la deviazione assoluta più significativa dalla media, non semplicemente la volatilità più alta.
1. Per ottenere la deviazione assoluta, sottraiamo ciascun valore dalla media e applichiamo la funzione abs().
2. Usando sort_values() possiamo ordinare i risultati dal più piccolo al più grande. Poi selezioniamo i 2 valori più grandi, ottobre e dicembre.
#take abs value move from the mean#we see Dec and Oct are the biggest abs moves
fin =abs(final - final.mean())print(fin.sort_values())
Oct_value = fin[10]
Dec_value = fin[12]print('Extreme Dec value:', Dec_value)print('Extreme Oct value:', Oct_value)Offerto da
date
90.044271110.10677130.12523750.15554020.21614610.36766170.41927160.42021880.57552140.852509101.080729121.231771Name: close,dtype: float64
Extreme Dec value:1.231770833333333
Extreme Oct value:1.080729166666667Offerto da
Simulazione
Ora che sappiamo qual è la domanda, dobbiamo costruire il nostro “modello nullo”.
Ci sono diverse opzioni:
Modelli parametrici. Se avessimo una buona idea della distribuzione dei dati o semplicemente facessimo ipotesi su di essa, potremmo usare tecniche “classiche” di test d’ipotesi, t-test, X², ANOVA a una via, ecc. Questi modelli possono essere restrittivi e una sorta di scatola nera se il ricercatore non comprende appieno le loro ipotesi.
Simulazione diretta. Potremmo fare ipotesi sul processo generativo dei dati e simulare direttamente. Per esempio, potremmo specificare un modello ARMA per i dati finanziari che stiamo modellando e costruirlo deliberatamente senza stagionalità. Potrebbe essere una scelta ragionevole per il nostro problema. Tuttavia, se conoscessimo il processo generativo dei dati dell’S&P500, saremmo già ricchi!
Simulazione tramite ricampionamento. È l’approccio che adotteremo. Campionando ripetutamente a caso dal dataset esistente e rimescolando le etichette, possiamo rendere l’effetto osservato ugualmente probabile tra tutte le etichette nei nostri dati (nel nostro caso le etichette sono le date), ottenendo così il dataset nullo desiderato.
Il campionamento è un argomento vasto, ma ci concentreremo su una tecnica in particolare, la permutazione, o shuffling.
Per ottenere il modello nullo desiderato, dobbiamo costruire un dataset che non presenti stagionalità. Se la nulla è vera, cioè non c’è stagionalità nei dati e l’effetto osservato è avvenuto semplicemente per caso, allora le etichette di ciascun mese (gen, feb, ecc.) sono prive di significato e quindi possiamo rimescolare ripetutamente i dati per costruire quella che la statistica classica chiamerebbe la “distribuzione campionaria della statistica di test sotto l’ipotesi nulla”. Questo ha l’effetto desiderato di rendere il fenomeno osservato (il valore estremo di dicembre) ugualmente probabile per tutti i mesi, esattamente ciò che richiede il nostro modello nullo.
Per dimostrare quanto siano potenti le tecniche di simulazione con l’hardware moderno, il codice in questo esempio permuterà in realtà i dati dei prezzi giornalieri, il che richiede molta più potenza di calcolo, ma si completa comunque in pochi secondi su una CPU moderna.
Nota: rimescolare le etichette giornaliere o mensili, nel nostro caso, produrrà comunque il dataset nullo desiderato.
Shuffle 'date' label to create null dataset
Offerto da
Una grande risorsa per imparare il campionamento è di Julian Simon.
Nota: Il modo in cui è costruito il nostro test equivale a un test a due code usando metodi “classici”, come il t-test di Welch o ANOVA, ecc., perché siamo interessati al valore più estremo, sopra o sotto la media.
Queste decisioni sono scelte di design, e abbiamo questa libertà perché il modello nullo è proprio questo: un Modello! Significa che possiamo specificarne i parametri come scegliamo; la chiave è essere davvero chiari su quale domanda stiamo cercando di rispondere.
Vogliamo quindi usare Python per simulare molti dati e creare un dataset nullo. Per farlo, simuleremo 1000 insiemi di 12 AMVR, permutando ogni volta le etichette 'date' per costruire la distribuzione campionaria. L’output di questo codice è incluso più avanti nella sezione sul p-hacking.
1. Per prima cosa definiamo alcuni contenitori per memorizzare i risultati: useremo un pd.DataFrame() e un semplice array [ ].
2. Poi definiamo un contatore a zero e avviamo un for-loop di 1000 iterazioni per creare i dati simulati.
3. La prima riga nel for-loop prende i rendimenti giornalieri originali (dall’inizio del tutorial) e usa il metodo sample() di pandas. Questo campiona casualmente dai rendimenti giornalieri un dato numero di volte. Nel nostro caso, abbiamo 8191 punti dati (252 giorni di trading per 32 anni). Eliminiamo anche l’indice con reset_index(), per poter aggiungere un nuovo indice. Serve perché l’indice data originale verrebbe rimescolato insieme ai dati, mentre vogliamo “rimescolare” solo i dati.
4. La riga successiva aggiunge il nuovo indice data riassegnando l’indice di daily_ret_shuffle con un pd.bdate_range di lunghezza pari ai dati originali. Usiamo bdate_range invece di date_range, perché vogliamo giorni lavorativi (5, lun-ven) e non 7 giorni a settimana.
5. Ora che abbiamo generato i nostri dati “rimescolati”, eseguiamo lo stesso data wrangling fatto all’inizio per costruire i valori AMVR. È ciò che fanno le 3 righe successive.
6. Una volta ottenuti i nostri valori AMVR simulati, dobbiamo aggiungerli al dataframe che memorizzerà tutti i 1000 run della simulazione. pd.concat prende il dataframe corrente e accoda ciascun dataframe successivo. Scegliamo axis 1 per avere colonne di dati.
7. Maximonth è una variabile che memorizza solo il valore più alto di ciascuna simulazione (lo useremo più avanti nell’analisi).
8. Poiché la nostra analisi richiede valori AMVR assoluti, le 3 righe successive prendono tutti i 1000 run AMVR e li appiattiscono in un unico grande array. Poi calcoliamo la media e sottraiamo i singoli valori da questa media e infine prendiamo il valore assoluto.
9. Facciamo lo stesso per i dati “solo massimi”. Nota: qui abbiamo usato una lista anziché un dataframe, il che ci permette di usare una list comprehension per calcolare l’abs(AMVR) per ciascuno dei nostri valori massimi.
Ora abbiamo tutti i dati necessari per completare la nostra analisi, sia le osservazioni originali sia il nuovo dataset simulato.
#as our Null is that no seasonality exists or alternatively that the month does not matter in terms of AMVR,#we can shuffle 'date' labels#for simplicity, we will shuffle the 'daily' return data, which has the same effect as shuffling 'month' labels#generate null data
new_df_sim = pd.DataFrame()
highest_only =[]
count=0
n=1000for i inrange(n):#sample same size as dataset, drop timestamp
daily_ret_shuffle = daily_ret.sample(8191).reset_index(drop=True)#add new timestamp to shuffled data
daily_ret_shuffle.index =(pd.bdate_range(start='1986-1-3',periods=8191))#then follow same data wrangling as before...
mnthly_annu = daily_ret_shuffle.resample('M').std()* np.sqrt(12)
ranked = mnthly_annu.groupby(mnthly_annu.index.year).rank()
sim_final = ranked.groupby(ranked.index.month).mean()#add each of 1000 sims into df
new_df_sim = pd.concat([new_df_sim,sim_final],axis=1)#also record just highest AMVR for each year (we will use this later for p-hacking explanation)
maxi_month =max(sim_final)
highest_only.append(maxi_month)#calculate absolute deviation in AMVR from the mean
all_months = new_df_sim.values.flatten()
mu_all_months = all_months.mean()
abs_all_months =abs(all_months-mu_all_months)#calculate absolute deviation in highest only AMVR from the mean
mu_highest = np.mean(highest_only)
abs_highest =[abs(x - mu_all_months)for x in highest_only]Offerto da
p-Hacking
Ecco la parte interessante. Abbiamo costruito un’ipotesi da testare, abbiamo generato dati simulati rimescolando le etichette ‘date’ dei dati, ora dobbiamo eseguire il nostro test d’ipotesi per trovare la probabilità di osservare un risultato significativo quanto quello di dicembre dato che l’ipotesi nulla (assenza di stagionalità) è vera.
Prima di eseguire il test, fissiamo le aspettative.
Qual è la probabilità di vedere almeno un risultato significativo dato un livello di significatività del 5%?
= 1-p(non significativo)
= 1-(1–0.05)¹²
= 0.46
quindi c’è il 46% di probabilità di vedere almeno un mese con un risultato significativo, dato che la nulla è vera.
Ora chiediamoci, per ciascun test individuale (confrontando l’AMVR assoluto di ciascuno dei 12 mesi con la media) quanti valori significativi dovremmo attenderci tra i nostri dati casuali, non stagionali?
12 x 0.05 = 0.6
Quindi con un livello di significatività 0,05, dovremmo attenderci un tasso di falsi positivi pari a 0,6. In altre parole, per ogni test (con i dati della nulla) confrontando l’AMVR di tutti i 12 mesi con la media, 0,6 mesi mostreranno un risultato significativo. (ovviamente non possiamo avere meno di 1 mese con risultato, ma con test ripetuti la media tenderà a questo numero).
Per tutta questa analisi abbiamo sottolineato l’importanza di essere molto chiari sulla domanda a cui vogliamo rispondere. Il problema con le aspettative che abbiamo appena calcolato è che abbiamo assunto di testare un risultato significativo su tutti i 12 mesi! Ecco perché la probabilità di vedere almeno un falso positivo è così alta, al 46%.
Questo è un esempio di bias da confronti multipli, in cui abbiamo ampliato lo spazio di ricerca e aumentato la probabilità di trovare un risultato significativo. È un problema perché potremmo abusare di questo fenomeno per selezionare i parametri del nostro modello che ci danno il p-value “desiderato”.
Questa è l’essenza del p-hacking
Per illustrare l’effetto del p-hacking e come ridurre la molteplicità, dobbiamo capire la differenza sottile ma significativa tra le seguenti 2 domande:
“Qual è la probabilità che dicembre appaia così estremo per caso?”
“Qual è la probabilità che un qualsiasi mese appaia così estremo per caso?”
La bellezza della simulazione sta nella sua semplicità. Il codice seguente è tutto ciò che serve per calcolare il p-value per rispondere alla prima domanda. Contiamo semplicemente quanti valori nel nostro dataset, usando tutte le 12000 deviazioni AMVR (12 mesi x 1000 prove), sono maggiori del valore osservato di dicembre. Otteniamo un p-value del 4,4%, vicino al nostro arbitrario cutoff del 5%, ma comunque significativo.
#count number of months in sim data where ave-vol-rank is >= Dec#Note: we are using Dec not Oct, as Dec has highest absolute deviation from the mean
count=0for i in abs_all_months:if i> Dec_value:
count+=1
ans = count/len(abs_all_months)print('p-value:', ans )Offerto da
Per rispondere alla seconda domanda ed evitare la molteplicità, invece di confrontare il nostro risultato con la distribuzione costruita con tutte le 12000 deviazioni AMVR, consideriamo solo il valore più alto da ciascuna delle 1000 prove dell’AMVR assoluto. Questo dà un p-value del 23%, decisamente non significativo!
#same again but just considering highest AMVR for each of 100 trials
count=0for i in abs_highest:if i> Dec_value:
count+=1
ans = count/len(abs_highest)print('p-value:', ans )Offerto da
Ora che abbiamo i risultati finali, tracciamo queste distribuzioni per mostrare visivamente l’effetto del p-hacking e gli esiti della nostra analisi:
1. Per prima cosa usiamo np.quantile() per trovare il nostro livello di significatività al 5%; lo tracceremo con axvline nei grafici inferiori.
2. Poi definiamo i 4 sottografici. Qui plt.subplots restituisce 'fig' = la nostra figura e 'ax1,ax2,ax3,ax4' = variabili che rappresentano ciascuno dei 4 sottografici.
3. Il grafico 1 mostra la prima colonna. Definiamo semplicemente un istogramma dei nostri dati 'abs_all_months' e scegliamo type='bar'. Per i grafici inferiori impostiamo cumulative='True', ottenendo così una cdf anziché una pdf. Bins definisce quante barre abbiamo per riempire ciascun intervallo di valori; 30 è un numero ragionevole per i nostri dati. Formattiamo poi i grafici con tecniche già viste, come axvline per tracciare i valori di significatività e di osservazione.
<matplotlib.legend.Legend at 0x280b4eb0b00>Offerto da
La colonna di sinistra contiene i dati che rispondono alla domanda 1 e quella di destra alla domanda 2. La riga superiore mostra le distribuzioni di probabilità e quella inferiore le CDF. La linea tratteggiata rossa è il livello di significatività al 5% che abbiamo deciso arbitrariamente. La linea blu è il valore AMVR estremo originale di dicembre pari a 1,23.
Il grafico a sinistra mostra che il valore di dicembre originale è significativo al 5%, ma di poco! Tuttavia, quando teniamo conto del bias da confronti multipli, nel grafico di destra la soglia di significatività si sposta da circa 1,2 (AMVR assoluto) a circa 1,6 (vedi la linea rossa).
Tenendo conto del bias da confronti multipli, il nostro valore di dicembre pari a 1,23 non è più significativo!
Considerando la domanda specifica a cui vogliamo rispondere ed evitando il bias da confronti multipli, abbiamo evitato di fare p-hacking sul nostro modello ed evitare di mostrare un risultato significativo quando non lo è.
Per esplorare ulteriormente il p-hacking e come possa essere abusato per raccontare una certa storia sui nostri dati, prova questa ottima app interattiva di FiveThirtyEight
Conclusioni
Abbiamo imparato che il test d’ipotesi non è il mostro spaventoso che pensavamo. Segui semplicemente i 3 passi sopra per costruire il tuo modello per qualsiasi tipo di dato o statistica di test.
Abbiamo mostrato che porre la domanda giusta è vitale per l’analisi scientifica. Una piccola variazione nella formulazione può portare a un modello molto diverso con risultati molto diversi.
Si spera sia ora evidente la potenza delle funzionalità Python più intermedie e la loro capacità di permetterti di condurre test statistici! Con poche righe di codice, abbiamo costruito e testato un fenomeno reale e tratto conclusioni azionabili.
Abbiamo discusso l’importanza di riconoscere e correggere il bias da confronti multipli ed evitare le trappole del p-hacking, mostrando come un risultato apparentemente significativo possa diventare non significativo.
Con sempre più “big data”, insieme alla pressione accademica a produrre ricerche con risultati “nuovi” o la pressione politica per mostrare un risultato “significativo”, la tentazione del p-hacking è in aumento. Imparando a riconoscere quando ne siamo colpevoli e a correggerci di conseguenza, possiamo diventare ricercatori migliori e produrre risultati scientifici più accurati e quindi azionabili!
Note dell’autore: I nostri risultati differiscono leggermente dalla ricerca originale di Winton, in parte per via di un dataset diverso (32 anni contro 87) e loro hanno ottobre come mese di interesse mentre noi abbiamo dicembre. Inoltre, hanno usato un metodo non divulgato per i loro “dati simulati” mentre noi abbiamo reso esplicita, attraverso il codice, la nostra metodologia di creazione di quei dati. Abbiamo fatto alcune assunzioni di modellazione lungo tutto il lavoro; anche queste sono state rese esplicite e si possono vedere nel codice. Queste scelte di design e di modellazione fanno parte del processo scientifico: purché siano esplicite, l’analisi ha valore.
Scopri i 15 migliori server MCP remoti che stanno trasformando lo sviluppo AI nel 2026. Scopri come migliorano automazione, ragionamento, sicurezza e velocità dei workflow.
%22%3E%3Cpath%20d%3D%22m-369%200h408v79h-408z%22%20fill%3D%22%237933ff%22%2F%3E%3Cpath%20d%3D%22m-148%200h306v45h-306z%22%20fill%3D%22%23974dff%22%2F%3E%3C%2Fg%3E%3C%2Fsvg%3E)
%22%2F%3E%3C%2FclipPath%3E%3Cg%20clip-path%3D%22url(%23a)%22%3E%3Cpath%20d%3D%22m76%200h484v36h-484z%22%20fill%3D%22%23974dff%22%2F%3E%3Cpath%20d%3D%22m0%200h384v24h-384z%22%20fill%3D%22%2306bdfc%22%20transform%3D%22matrix(1%200%200%20-1%200%2060)%22%2F%3E%3C%2Fg%3E%3C%2Fsvg%3E)
Abid Ali Awan · 10 min
%22%3E%0A%20%20%20%20%20%20%20%20%3Crect%20fill%3D%22%2305192D%22%20height%3D%22180%22%20width%3D%22336%22%20%2F%3E%0A%20%20%20%20%20%20%20%20%3Crect%0A%20%20%20%20%20%20%20%20%20%20fill%3D%22%23B08EFF%22%0A%20%20%20%20%20%20%20%20%20%20height%3D%22106.6%22%0A%20%20%20%20%20%20%20%20%20%20transform%3D%22rotate(-30%20-208.588%20220.18)%22%0A%20%20%20%20%20%20%20%20%20%20width%3D%22456.125%22%0A%20%20%20%20%20%20%20%20%20%20x%3D%22-208.588%22%0A%20%20%20%20%20%20%20%20%20%20y%3D%22220.18%22%0A%20%20%20%20%20%20%20%20%2F%3E%0A%20%20%20%20%20%20%20%20%3Crect%0A%20%20%20%20%20%20%20%20%20%20fill%3D%22%237933FF%22%0A%20%20%20%20%20%20%20%20%20%20height%3D%2235.875%22%0A%20%20%20%20%20%20%20%20%20%20transform%3D%22rotate(-30%20-208.588%20220.18)%22%0A%20%20%20%20%20%20%20%20%20%20width%3D%22347.576%22%0A%20%20%20%20%20%20%20%20%20%20x%3D%22-208.588%22%0A%20%20%20%20%20%20%20%20%20%20y%3D%22220.18%22%0A%20%20%20%20%20%20%20%20%2F%3E%0A%20%20%20%20%20%20%20%20%3Crect%0A%20%20%20%20%20%20%20%20%20%20fill%3D%22%2372E5FE%22%0A%20%20%20%20%20%20%20%20%20%20height%3D%2235.875%22%0A%20%20%20%20%20%20%20%20%20%20transform%3D%22rotate(-30%20-226.525%20189.111)%22%0A%20%20%20%20%20%20%20%20%20%20width%3D%22347.576%22%0A%20%20%20%20%20%20%20%20%20%20x%3D%22-226.525%22%0A%20%20%20%20%20%20%20%20%20%20y%3D%22189.111%22%0A%20%20%20%20%20%20%20%20%2F%3E%0A%20%20%20%20%20%20%20%20%3Crect%0A%20%20%20%20%20%20%20%20%20%20fill%3D%22%235EB1FF%22%0A%20%20%20%20%20%20%20%20%20%20height%3D%22106.6%22%0A%20%20%20%20%20%20%20%20%20%20transform%3D%22rotate(-30%2085.5879%20256.055)%22%0A%20%20%20%20%20%20%20%20%20%20width%3D%22456.125%22%0A%20%20%20%20%20%20%20%20%20%20x%3D%2285.5879%22%0A%20%20%20%20%20%20%20%20%20%20y%3D%22256.055%22%0A%20%20%20%20%20%20%20%20%2F%3E%0A%20%20%20%20%20%20%20%20%3Crect%0A%20%20%20%20%20%20%20%20%20%20fill%3D%22%23B08EFF%22%0A%20%20%20%20%20%20%20%20%20%20height%3D%2235.875%22%0A%20%20%20%20%20%20%20%20%20%20transform%3D%22rotate(-30%2085.5879%20256.055)%22%0A%20%20%20%20%20%20%20%20%20%20width%3D%22184.5%22%0A%20%20%20%20%20%20%20%20%20%20x%3D%2285.5879%22%0A%20%20%20%20%20%20%20%20%20%20y%3D%22256.055%22%0A%20%20%20%20%20%20%20%20%2F%3E%0A%20%20%20%20%20%20%20%20%3Crect%0A%20%20%20%20%20%20%20%20%20%20fill%3D%22%2372E5FE%22%0A%20%20%20%20%20%20%20%20%20%20height%3D%2235.875%22%0A%20%20%20%20%20%20%20%20%20%20transform%3D%22rotate(-30%20227.432%20132.736)%22%0A%20%20%20%20%20%20%20%20%20%20width%3D%22271.625%22%0A%20%20%20%20%20%20%20%20%20%20x%3D%22227.432%22%0A%20%20%20%20%20%20%20%20%20%20y%3D%22132.736%22%0A%20%20%20%20%20%20%20%20%2F%3E%0A%20%20%20%20%20%20%3C%2Fg%3E%0A%20%20%20%20%20%20%3Cdefs%3E%0A%20%20%20%20%20%20%20%20%3CclipPath%20id%3D%22clip%22%3E%0A%20%20%20%20%20%20%20%20%20%20%3Crect%20fill%3D%22white%22%20height%3D%22180%22%20width%3D%22336%22%20%2F%3E%0A%20%20%20%20%20%20%20%20%3C%2FclipPath%3E%0A%20%20%20%20%20%20%3C%2Fdefs%3E%0A%20%20%20%20%3C%2Fsvg%3E)
Tim Lu · 12 min
"Negli ultimi 32 anni, ottobre è stato in media il mese più volatile per l’S&P500 e dicembre il meno volatile".
In questo tutorial useremo Python per eseguire un’analisi completa e un test di questo fenomeno, per capire se è statisticamente significativo oppure no.
Useremo le funzionalità per serie temporali di pandas, come Resample, per trasformare i prezzi azionari grezzi in un formato che dia vita ai dati.
Scopriremo poi come condurre test d’ipotesi usando la simulazione in Python invece di formule macchinose.
Infine parleremo di un problema importante nell’analisi statistica, il bias da confronti multipli, vedremo come gestirlo e mostreremo visivamente con matplotlib l’effetto del “p-Hacking”.
Il nostro obiettivo:
I nostri dati:
Useremo dati giornalieri dell’S&P500 per questa analisi; in particolare, useremo i prezzi di chiusura giornalieri grezzi dal 1986 al 2018 (sorprendentemente difficili da trovare, quindi li ho resi pubblicamente disponibili).
L’ispirazione per questo post viene da Winton, che riprodurremo qui, seppur con 32 anni di dati contro i loro 87.
Lavorare i dati con Pandas
Per rispondere alla domanda se l’elevata volatilità osservata in alcuni mesi sia davvero significativa e quindi probabile che continui, dobbiamo trasformare i nostri 32 anni di dati sui prezzi in un formato che evidenzi il fenomeno che stiamo indagando.
Il formato scelto sarà il ranking medio mensile della volatilità (AMVR).
Il codice seguente mostra come portiamo i dati grezzi dei prezzi in questo formato. Iniziamo!
Per prima cosa le importazioni standard. (matplotlib.patches ci permette di controllare lo stile delle singole barre in un istogramma)
Un trucco utile per far sì che i grafici occupino tutta la larghezza nei notebook Jupyter:
Questa è opzionale: provala e guarda cosa succede!
Importiamo i dati usando il metodo read_csv, che prende un percorso, nel nostro caso un URL ai dati. Indichiamo anche di usare "date" come indice e di analizzare automaticamente le date dal testo fornito, al meglio possibile.
Questo ci dà i prezzi di chiusura grezzi e non rettificati dell’S&P500 (SPX). Ora dobbiamo convertire i prezzi grezzi in rendimenti giornalieri in percentuale. Per farlo, abbiamo due opzioni:
Per il nostro scopo useremo il metodo 2, dato che pandas può gestire istantaneamente il calcolo su un dataset di questa dimensione (oltre 8000 valori). Puliremo anche i dati eliminando il primo valore che diventa 'NaN' perché non c’è variazione di prezzo rispetto al giorno precedente. pct_change() accetta un parametro opzionale 'periods' che modifica lo shift temporale. Lo lasciamo all’impostazione predefinita pari a uno.
Il passo successivo è prendere queste variazioni giornaliere in percentuale e trasformarle in “volatilità mensile annualizzata”. La prima riga di codice sotto mostra che possiamo effettuare questa trasformazione in una sola riga. Questa è la potenza di Pandas! Ma scomponiamo ogni passaggio in sequenza.
1. Per ottenere mnthly_annu, usiamo innanzitutto il metodo 'resample' sui nostri rendimenti giornalieri. Resample ci permette di cambiare la frequenza dei periodi. Accetta una 'frequency offset string' o, in parole semplici, una lettera che corrisponde alla nuova frequenza desiderata. Includono:
Poiché vogliamo passare da rendimenti giornalieri a mensili, passiamo 'M' come parametro.
2. La seconda cosa è decidere come arrivare al nuovo valore mensile; ad esempio, potremmo sommare i valori giornalieri, moltiplicarli, ecc. Per la nostra analisi vogliamo una misura di volatilità e la deviazione standard funziona bene, quindi aggiungeremo std() ai nostri dati ricampionati per ottenere la volatilità mensile.
3. L’ultimo passo è annualizzare questo valore. Lo facciamo semplicemente moltiplicando per la radice quadrata di 12, con 12 che è il numero di periodi (mesi) in un anno. Questo ci dà i valori di volatilità mensile annualizzata che ci servono.
La buona pratica nell’analisi dei dati è visualizzarli man mano che si procede. Diamo quindi un’occhiata ai nostri dati di volatilità annualizzata.
Il grafico sotto mostra chiaramente i principali eventi di mercato, come il Black Monday e la Crisi Finanziaria del 2008. Il metodo matplotlib 'axvspan' ci consente di aggiungere linee verticali (axhspan è il corrispondente per le linee orizzontali); accetta 'xmin' e 'xmax' come parametri, che specificano l’ampiezza del rettangolo creato. Poiché il nostro asse x è l’indice datetime, possiamo semplicemente passare gli anni che vogliamo evidenziare. Il parametro alpha regola la trasparenza del blocco di colore per poter vedere il grafico sottostante.
L’attributo mpatches ci permette di creare una legenda personalizzata. Lo facciamo definendo prima la variabile 'labs' con colore, alpha e il testo desiderato, quindi passandola al parametro 'handles' in plt.legend per renderizzare la legenda.
Abbiamo quindi visto uno dei metodi potenti di pandas, resample. Ora usiamone un altro, 'groupby'. Dobbiamo passare dai valori di volatilità mensile annualizzata al nostro AMVR desiderato. Anche qui possiamo riuscirci in poche righe di codice.
1. Per prima cosa applichiamo 'groupby' ai valori di mnthly_annu. Groupby richiede un parametro che specifica come raggruppare i dati. Può essere una funzione o, come nel nostro caso, una Series. Passiamo 'mnthly_annu.index.year', che è semplicemente l’attributo anno dell’indice datetime di mnthly_annu. Questo raggrupperà i valori di volatilità mensile per ognuno dei 32 anni del nostro set di dati.
2. Successivamente, applichiamo il metodo 'rank', che ordina i dati in ordine crescente.
3. Infine, ripetiamo lo stesso e facciamo la media su tutti gli anni per ciascun mese. Questo ci dà i nostri valori AMVR finali!
Questo ci dà il nostro ranking medio mensile della volatilità. Numericamente vediamo che il mese 10 (ottobre) è il più alto e il 12 (dicembre) è il più basso.
Scegliere la giusta visualizzazione è importante perché vogliamo raccontare la storia dei dati con i grafici più efficaci. Per la nostra analisi vogliamo mostrare chiaramente la volatilità più alta e quella più bassa. Useremo quindi il grafico a barre di matplotlib, annotando e colorando di conseguenza per dare il massimo impatto.
1. Indicizzando b_plot (b_plot[9]) con la barra desiderata, possiamo impostarne il colore per evidenziare i valori massimo e minimo.
2. Per aggiungere i valori AMVR alla rispettiva barra, iteriamo su ciascun valore AMVR (variabile 'final') e arrotondiamo a 2 decimali. 'i' e 'v' rappresentano 'indice' e 'valore' per la variabile 'final', quindi 'i' andrà da 1 a 12 e 'value' dà il valore AMVR arrotondato a 2 decimali. Usiamo queste variabili in plt.text(). Il primo parametro di plt.text() è la posizione sull’asse x per l’etichetta, quindi passiamo 'i' per ogni barra e lo spostiamo di 0.8 per centrarla. Il secondo parametro è la stringa del valore di 'final'.
3. Per mostrare il valore medio, usiamo axhline e il parametro 'label' per aggiungere una legenda nell’angolo in alto a destra. Passiamo semplicemente il valore sull’asse y, cioè la media, e poi aggiungiamo lo stile.
Questo offre la visuale migliore del fenomeno che stiamo indagando. Si vede chiaramente che ottobre è stato il mese più volatile e dicembre il meno. Importante, si nota anche che dicembre è il valore più “estremo” in termini assoluti. Questo sarà importante per l’analisi nella sezione successiva.
Questi sono i nostri dati; passiamo al test d’ipotesi…
Test d’ipotesi: qual è la domanda?
Il test d’ipotesi è una delle tecniche fondamentali della data science, eppure è tra le più intimidatorie e fraintese. La base di questa paura sta in come viene insegnato in Statistica 101, dove ci viene detto di:
"eseguire un t-test; è monodirezionale o bidirezionale? Scegli un test statistico adatto come il t-test di Welch, calcola i gradi di libertà, calcola il punteggio t, cerca il valore critico in tabella, confronta il valore critico con la statistica t…"
Comprensibilmente, questo porta a confusione su quale test eseguire e come eseguirlo. Tuttavia, tutte queste tecniche statistiche classiche per eseguire un test d’ipotesi sono state sviluppate quando avevamo pochissima potenza di calcolo ed erano semplicemente soluzioni analitiche chiuse per calcolare un p-value, tutto qui! Ma con l’ulteriore complicazione di dover scegliere la formula giusta per la situazione data, a causa delle loro ipotesi restrittive e talvolta opache.
Ma rallegriamoci!
C’è un modo migliore. La simulazione.
Per capire come la simulazione possa aiutarci, ricordiamo cos’è un test d’ipotesi:
Vogliamo testare “se l’effetto osservato nei nostri dati sia reale o se possa avvenire semplicemente per caso” e per eseguire questo test facciamo quanto segue:
Esiste un solo test! Questo ottimo blog di Allen Downey approfondisce il test d’ipotesi
Il vero potere della simulazione è che dobbiamo rendere esplicite, attraverso il codice, le ipotesi del nostro modello. Mentre le tecniche classiche possono essere una “scatola nera” quanto alle loro ipotesi.
Esempio sotto: Il grafico a sinistra mostra i dati reali e l’effetto osservato con una certa probabilità (verde). Quello a destra è il nostro dato nullo simulato con la registrazione di quando l’effetto osservato è apparso per caso (rosso). Questa è la base del test d’ipotesi: qual è la probabilità di vedere l’effetto osservato sui dati della nulla.
La parte più critica del test d’ipotesi è chiarire che domanda stiamo cercando di rispondere. Nel nostro caso chiediamo:
“Il valore più estremo potrebbe verificarsi per caso?”
Il valore più estremo lo definiamo come la massima deviazione assoluta dell’AMVR dalla media. Questa domanda costituisce la nostra ipotesi nulla.
Nei nostri dati, il valore più estremo è quello di dicembre (1,23) non quello di ottobre (1,08), perché stiamo guardando la deviazione assoluta più significativa dalla media, non semplicemente la volatilità più alta.
1. Per ottenere la deviazione assoluta, sottraiamo ciascun valore dalla media e applichiamo la funzione abs().
2. Usando sort_values() possiamo ordinare i risultati dal più piccolo al più grande. Poi selezioniamo i 2 valori più grandi, ottobre e dicembre.
Simulazione
Ora che sappiamo qual è la domanda, dobbiamo costruire il nostro “modello nullo”.
Ci sono diverse opzioni:
Il campionamento è un argomento vasto, ma ci concentreremo su una tecnica in particolare, la permutazione, o shuffling.
Per ottenere il modello nullo desiderato, dobbiamo costruire un dataset che non presenti stagionalità. Se la nulla è vera, cioè non c’è stagionalità nei dati e l’effetto osservato è avvenuto semplicemente per caso, allora le etichette di ciascun mese (gen, feb, ecc.) sono prive di significato e quindi possiamo rimescolare ripetutamente i dati per costruire quella che la statistica classica chiamerebbe la “distribuzione campionaria della statistica di test sotto l’ipotesi nulla”. Questo ha l’effetto desiderato di rendere il fenomeno osservato (il valore estremo di dicembre) ugualmente probabile per tutti i mesi, esattamente ciò che richiede il nostro modello nullo.
Per dimostrare quanto siano potenti le tecniche di simulazione con l’hardware moderno, il codice in questo esempio permuterà in realtà i dati dei prezzi giornalieri, il che richiede molta più potenza di calcolo, ma si completa comunque in pochi secondi su una CPU moderna.
Nota: rimescolare le etichette giornaliere o mensili, nel nostro caso, produrrà comunque il dataset nullo desiderato.
Una grande risorsa per imparare il campionamento è di Julian Simon.
Nota: Il modo in cui è costruito il nostro test equivale a un test a due code usando metodi “classici”, come il t-test di Welch o ANOVA, ecc., perché siamo interessati al valore più estremo, sopra o sotto la media.
Queste decisioni sono scelte di design, e abbiamo questa libertà perché il modello nullo è proprio questo: un Modello! Significa che possiamo specificarne i parametri come scegliamo; la chiave è essere davvero chiari su quale domanda stiamo cercando di rispondere.
Vogliamo quindi usare Python per simulare molti dati e creare un dataset nullo. Per farlo, simuleremo 1000 insiemi di 12 AMVR, permutando ogni volta le etichette 'date' per costruire la distribuzione campionaria. L’output di questo codice è incluso più avanti nella sezione sul p-hacking.
1. Per prima cosa definiamo alcuni contenitori per memorizzare i risultati: useremo un pd.DataFrame() e un semplice array [ ].
2. Poi definiamo un contatore a zero e avviamo un for-loop di 1000 iterazioni per creare i dati simulati.
3. La prima riga nel for-loop prende i rendimenti giornalieri originali (dall’inizio del tutorial) e usa il metodo sample() di pandas. Questo campiona casualmente dai rendimenti giornalieri un dato numero di volte. Nel nostro caso, abbiamo 8191 punti dati (252 giorni di trading per 32 anni). Eliminiamo anche l’indice con reset_index(), per poter aggiungere un nuovo indice. Serve perché l’indice data originale verrebbe rimescolato insieme ai dati, mentre vogliamo “rimescolare” solo i dati.
4. La riga successiva aggiunge il nuovo indice data riassegnando l’indice di daily_ret_shuffle con un pd.bdate_range di lunghezza pari ai dati originali. Usiamo bdate_range invece di date_range, perché vogliamo giorni lavorativi (5, lun-ven) e non 7 giorni a settimana.
5. Ora che abbiamo generato i nostri dati “rimescolati”, eseguiamo lo stesso data wrangling fatto all’inizio per costruire i valori AMVR. È ciò che fanno le 3 righe successive.
6. Una volta ottenuti i nostri valori AMVR simulati, dobbiamo aggiungerli al dataframe che memorizzerà tutti i 1000 run della simulazione. pd.concat prende il dataframe corrente e accoda ciascun dataframe successivo. Scegliamo axis 1 per avere colonne di dati.
7. Maximonth è una variabile che memorizza solo il valore più alto di ciascuna simulazione (lo useremo più avanti nell’analisi).
8. Poiché la nostra analisi richiede valori AMVR assoluti, le 3 righe successive prendono tutti i 1000 run AMVR e li appiattiscono in un unico grande array. Poi calcoliamo la media e sottraiamo i singoli valori da questa media e infine prendiamo il valore assoluto.
9. Facciamo lo stesso per i dati “solo massimi”. Nota: qui abbiamo usato una lista anziché un dataframe, il che ci permette di usare una list comprehension per calcolare l’abs(AMVR) per ciascuno dei nostri valori massimi.
Ora abbiamo tutti i dati necessari per completare la nostra analisi, sia le osservazioni originali sia il nuovo dataset simulato.
p-Hacking
Ecco la parte interessante. Abbiamo costruito un’ipotesi da testare, abbiamo generato dati simulati rimescolando le etichette ‘date’ dei dati, ora dobbiamo eseguire il nostro test d’ipotesi per trovare la probabilità di osservare un risultato significativo quanto quello di dicembre dato che l’ipotesi nulla (assenza di stagionalità) è vera.
Prima di eseguire il test, fissiamo le aspettative.
Qual è la probabilità di vedere almeno un risultato significativo dato un livello di significatività del 5%?
= 1-p(non significativo)
= 1-(1–0.05)¹²
= 0.46
quindi c’è il 46% di probabilità di vedere almeno un mese con un risultato significativo, dato che la nulla è vera.
Ora chiediamoci, per ciascun test individuale (confrontando l’AMVR assoluto di ciascuno dei 12 mesi con la media) quanti valori significativi dovremmo attenderci tra i nostri dati casuali, non stagionali?
12 x 0.05 = 0.6
Quindi con un livello di significatività 0,05, dovremmo attenderci un tasso di falsi positivi pari a 0,6. In altre parole, per ogni test (con i dati della nulla) confrontando l’AMVR di tutti i 12 mesi con la media, 0,6 mesi mostreranno un risultato significativo. (ovviamente non possiamo avere meno di 1 mese con risultato, ma con test ripetuti la media tenderà a questo numero).
Per tutta questa analisi abbiamo sottolineato l’importanza di essere molto chiari sulla domanda a cui vogliamo rispondere. Il problema con le aspettative che abbiamo appena calcolato è che abbiamo assunto di testare un risultato significativo su tutti i 12 mesi! Ecco perché la probabilità di vedere almeno un falso positivo è così alta, al 46%.
Questo è un esempio di bias da confronti multipli, in cui abbiamo ampliato lo spazio di ricerca e aumentato la probabilità di trovare un risultato significativo. È un problema perché potremmo abusare di questo fenomeno per selezionare i parametri del nostro modello che ci danno il p-value “desiderato”.
Per illustrare l’effetto del p-hacking e come ridurre la molteplicità, dobbiamo capire la differenza sottile ma significativa tra le seguenti 2 domande:
La bellezza della simulazione sta nella sua semplicità. Il codice seguente è tutto ciò che serve per calcolare il p-value per rispondere alla prima domanda. Contiamo semplicemente quanti valori nel nostro dataset, usando tutte le 12000 deviazioni AMVR (12 mesi x 1000 prove), sono maggiori del valore osservato di dicembre. Otteniamo un p-value del 4,4%, vicino al nostro arbitrario cutoff del 5%, ma comunque significativo.
Per rispondere alla seconda domanda ed evitare la molteplicità, invece di confrontare il nostro risultato con la distribuzione costruita con tutte le 12000 deviazioni AMVR, consideriamo solo il valore più alto da ciascuna delle 1000 prove dell’AMVR assoluto. Questo dà un p-value del 23%, decisamente non significativo!
Ora che abbiamo i risultati finali, tracciamo queste distribuzioni per mostrare visivamente l’effetto del p-hacking e gli esiti della nostra analisi:
1. Per prima cosa usiamo np.quantile() per trovare il nostro livello di significatività al 5%; lo tracceremo con axvline nei grafici inferiori.
2. Poi definiamo i 4 sottografici. Qui plt.subplots restituisce 'fig' = la nostra figura e 'ax1,ax2,ax3,ax4' = variabili che rappresentano ciascuno dei 4 sottografici.
3. Il grafico 1 mostra la prima colonna. Definiamo semplicemente un istogramma dei nostri dati 'abs_all_months' e scegliamo type='bar'. Per i grafici inferiori impostiamo cumulative='True', ottenendo così una cdf anziché una pdf. Bins definisce quante barre abbiamo per riempire ciascun intervallo di valori; 30 è un numero ragionevole per i nostri dati. Formattiamo poi i grafici con tecniche già viste, come axvline per tracciare i valori di significatività e di osservazione.
La colonna di sinistra contiene i dati che rispondono alla domanda 1 e quella di destra alla domanda 2. La riga superiore mostra le distribuzioni di probabilità e quella inferiore le CDF. La linea tratteggiata rossa è il livello di significatività al 5% che abbiamo deciso arbitrariamente. La linea blu è il valore AMVR estremo originale di dicembre pari a 1,23.
Il grafico a sinistra mostra che il valore di dicembre originale è significativo al 5%, ma di poco! Tuttavia, quando teniamo conto del bias da confronti multipli, nel grafico di destra la soglia di significatività si sposta da circa 1,2 (AMVR assoluto) a circa 1,6 (vedi la linea rossa).
Considerando la domanda specifica a cui vogliamo rispondere ed evitando il bias da confronti multipli, abbiamo evitato di fare p-hacking sul nostro modello ed evitare di mostrare un risultato significativo quando non lo è.
Per esplorare ulteriormente il p-hacking e come possa essere abusato per raccontare una certa storia sui nostri dati, prova questa ottima app interattiva di FiveThirtyEight
Conclusioni
Note dell’autore: I nostri risultati differiscono leggermente dalla ricerca originale di Winton, in parte per via di un dataset diverso (32 anni contro 87) e loro hanno ottobre come mese di interesse mentre noi abbiamo dicembre. Inoltre, hanno usato un metodo non divulgato per i loro “dati simulati” mentre noi abbiamo reso esplicita, attraverso il codice, la nostra metodologia di creazione di quei dati. Abbiamo fatto alcune assunzioni di modellazione lungo tutto il lavoro; anche queste sono state rese esplicite e si possono vedere nel codice. Queste scelte di design e di modellazione fanno parte del processo scientifico: purché siano esplicite, l’analisi ha valore.
Se vuoi saperne di più sulla finanza in Python, segui il corso di DataCamp Intro to Python for Finance e dai un’occhiata al Tutorial Python per la finanza: trading algoritmico.
Seguimi su twitter.com/pdquant per altro!