Leer tijdreeksgegevens manipuleren met pandas en significantietoetsen uitvoeren via simulatie met Python, om de volatiliteit van de aandelenmarkt te analyseren.
Oktober is historisch gezien de meest volatiele maand voor aandelen, maar is dit een blijvend signaal of slechts ruis in de data?
"In de afgelopen 32 jaar was oktober gemiddeld de meest volatiele maand voor de S&P500 en december de minst volatiele".
In deze tutorial gebruiken we Python om een volledige analyse en toetsing van dit fenomeen uit te voeren en te bepalen of het statistisch significant is of niet.
We gebruiken tijdreeksfunctionaliteit van pandas, zoals Resample, om ruwe aandelenkoersdata te bewerken tot een formaat dat de data tot leven brengt.
Vervolgens ontdekken we hoe je hypothesetoetsing uitvoert met simulatie in Python in plaats van omslachtige formules.
Tot slot bespreken we een groot probleem in statistische analyse, Multiple Comparison Bias, leren we hoe je daarmee omgaat en laten we via matplotlib visueel het effect van 'p-hacking' zien.
Ons doel:
Demonstreren hoe je Pandas gebruikt om tijdreeksen te analyseren
Begrijpen hoe je een hypothesetoets opstelt
Simulatie met Python gebruiken om hypothesetoetsing uit te voeren
Het belang tonen van corrigeren voor multiple comparison bias
Onze data:
We gebruiken dagelijkse S&P500-data voor deze analyse; specifiek gebruiken we de ruwe dagelijkse slotkoersen van 1986 tot 2018 (die verrassend moeilijk te vinden zijn, dus ik heb ze openbaar beschikbaar gemaakt).
De inspiratie voor dit artikel kwam van Winton, wat we hier zullen reproduceren, zij het met 32 jaar data tegenover hun 87 jaar.
Aan de slag met Pandas
Om de vraag te beantwoorden of de extreme volatiliteit in bepaalde maanden echt significant is en dus waarschijnlijk zal aanhouden, moeten we onze 32 jaar aan prijsdata omzetten naar een formaat dat het fenomeen dat we onderzoeken zichtbaar maakt.
Ons voorkeursformaat is de gemiddelde maandelijkse volatiliteitsrangschikking (AMVR).
De volgende code laat zien hoe we onze ruwe prijsdata in dit formaat krijgen. Aan de slag!
Eerst de standaard imports. (matplotlib.patches geeft ons controle over de styling van individuele balken in een histogram)
#standard importsimport numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
%matplotlib inline
Mogelijk gemaakt door
Een handige truc om grafieken op volledige breedte weer te geven in Jupyter-notebooks:
#resize charts to fit screen if using Jupyter Notebook
plt.rcParams['figure.figsize']=[15,5]Mogelijk gemaakt door
Deze is optioneel - probeer het eens en kijk wat er gebeurt!
Importeer onze data met de methode read_csv, die een pad krijgt, in ons geval een URL naar de data. We geven ook aan dat 'date' als index moet worden gebruikt en dat de datums automatisch zo goed mogelijk uit de tekst moeten worden geparseerd.
#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()Mogelijk gemaakt door
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
Dit levert de onbewerkte, niet-aangepaste slotkoersen van de S&P500 (SPX) op. Nu moeten we onze ruwe prijzen omzetten in dagelijkse %‑rendementen. Daarvoor hebben we twee opties:
We kunnen de natuurlijke logaritme van de prijzen nemen. Dit geeft een benadering van de werkelijke dagelijkse rendementen.
We kunnen de pandas‑methode 'pct_change()' gebruiken om de dagelijkse procentuele verandering direct te berekenen.
Voor ons doel gebruiken we methode 2, omdat pandas de berekening direct aankan op een dataset van deze omvang (meer dan 8000 waarden). We schonen de data ook op door de eerste waarde te verwijderen, die 'NaN' wordt omdat er geen prijsverandering is ten opzichte van de dag ervoor. pct_change() neemt een optionele parameter 'periods' die de periodeshift verandert. We laten die op de standaardwaarde één staan.
#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()Mogelijk gemaakt door
date
1986-01-030.0061551986-01-06-0.0010911986-01-070.0149541986-01-08-0.0272681986-01-09-0.008944Name: close,dtype: float64
Mogelijk gemaakt door
De volgende stap is om deze dagelijkse %‑veranderingen te transformeren naar 'maandelijkse geannualiseerde volatiliteit'. De eerste regel code hieronder laat zien dat we deze transformatie in één regel kunnen uitvoeren. Dit is de kracht van pandas! Maar laten we elke stap op volgorde uitpakken.
1. Om mnthly_annu te krijgen, gebruiken we eerst de methode 'resample' op onze dagelijkse rendementen. Resample laat ons de frequentie van de dataperiodes veranderen. Het neemt een 'frequency offset string' of, simpel gezegd, een letter die overeenkomt met de gewenste nieuwe frequentie. Deze omvatten:
B business day‑frequentie
C custom business day‑frequentie
D kalenderdag‑frequentie
W wekelijkse frequentie
M maandultimo‑frequentie
Aangezien we willen resamplen van dagelijkse rendementen naar maandelijks, geven we 'M' door als parameter.
2. Ten tweede moeten we beslissen hoe we tot deze nieuwe maandwaarde willen komen; we zouden de dagelijkse waarden kunnen sommeren, vermenigvuldigen, enzovoort. Voor onze analyse willen we een maat voor volatiliteit en standaarddeviatie werkt goed, dus voegen we std() toe aan onze geresamplede data om maandelijkse volatiliteit te krijgen.
3. De laatste stap is om dit cijfer te annualiseren. Dat doen we eenvoudig door te vermenigvuldigen met de wortel uit 12, waarbij 12 het aantal periodes (maanden) per jaar is. Dit levert de geannualiseerde maandvolatiliteiten op die we nodig hebben.
Best practice in data‑analyse is om de data te visualiseren terwijl we door de analyse heen gaan. Laten we dus naar onze geannualiseerde volatiliteitsdata kijken.
De onderstaande plot toont duidelijk grote marktgebeurtenissen, zoals Black Monday en de financiële crisis van 2008. De matplotlib‑methode 'axvspan' laat ons verticale lijnen toevoegen (axhspan is de overeenkomstige methode voor horizontale lijnen). Hij neemt 'xmin' en 'xmax' als parameters, die de breedte van de gecreëerde rechthoek specificeren; omdat onze x‑as de datetime‑index is, kunnen we eenvoudig de jaren doorgeven die we willen markeren. De alpha‑parameter past de transparantie van de blokkleur aan zodat we de grafiek eronder nog kunnen zien.
Het attribuut mpatches laat ons een aangepaste legenda maken. Dat doen we door eerst de variabele 'labs' te definiëren met kleur, alpha en de gewenste tekst, en die vervolgens door te geven aan de 'handles'‑parameter binnen plt.legend om de legenda te renderen.
#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])Mogelijk gemaakt door
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>Mogelijk gemaakt door
We hebben nu een van de krachtige methodes in pandas gezien: resample. Laten we nu een andere gebruiken, 'groupby'. We moeten van onze maandelijkse geannualiseerde volatiliteitswaarden naar onze gewenste AMVR‑maatstaf. Ook dit kunnen we in slechts een paar regels code bereiken.
1. Eerst passen we de 'groupby'‑methode toe op onze mnthly_annu‑waarden. Groupby vereist een parameter die specificeert hoe de data gegroepeerd moet worden. Dit kan een functie zijn of, zoals in ons geval, een Series. We geven 'mnthly_annu.index.year' door; dit is simpelweg het jaarattribuut op de datetime‑index van de mnthly_annu‑data. Dit groepeert onze maandvolatiliteitswaarden voor elk van de 32 jaar in onze dataset.
2. Vervolgens passen we de methode 'rank' toe, die de data in oplopende volgorde rangschikt.
3. Tot slot doen we hetzelfde nog eens en middelen we over alle jaren per maand. Dit levert onze definitieve AMVR‑waarden op!
#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()Mogelijk gemaakt door
count 12.000000
mean 6.450521
std 0.627458
min 5.21875025%6.03101350%6.49100475%6.704545
max 7.531250Name: close,dtype: float64
Mogelijk gemaakt door
Dit levert onze definitieve Gemiddelde Maandelijkse Volatiliteitsrangschikking op. Numeriek zien we dat maand 10 (oktober) het hoogst is en 12 (december) het laagst.
date
16.81818226.66666736.57575847.30303056.60606166.03030376.03125085.87500096.406250107.531250116.343750125.218750Name: close,dtype: float64
Mogelijk gemaakt door
De juiste visualisatie kiezen voor onze data is belangrijk, omdat we het verhaal van de data willen vertellen met de meest impactvolle grafieken en plots. Voor onze analyse willen we duidelijk de hoogste en de laagste volatiliteit tonen. Daarvoor gebruiken we de staafdiagram van matplotlib en annoteren en kleuren we dienovereenkomstig voor het meeste effect.
1. Door b_plot (b_plot[9]) te indexeren met de gewenste balk kunnen we de kleur instellen om de hoogste en laagste waarden te markeren.
2. Om de AMVR‑waarden aan hun respectieve balk toe te voegen, lopen we door elke AMVR‑waarde ('final'‑variabele) en ronden we het resultaat af op 2 decimalen. 'i' en 'v' staan voor de 'index' en 'waarde' van de 'final'‑variabele, dus 'i' loopt van 1–12 en 'value' geeft de AMVR‑waarde afgerond op 2 decimalen. We gebruiken deze variabelen binnen plt.text(). De eerste parameter van plt.txt() is de x‑aslocatie voor het label, dus we geven 'i' door voor elke balk en verschuiven met 0,8 om te centreren. De tweede parameter is de stringversie van de 'final'‑waarde.
3. Om de gemiddelde waarde te tonen gebruiken we axhline en gebruiken we de 'label'‑waarde om rechtsboven een legenda toe te voegen. Geef simpelweg de y‑aswaarde door, oftewel het gemiddelde, en voeg vervolgens styling toe.
Dit levert de beste visual tot nu toe van het fenomeen dat we onderzoeken. We zien duidelijk dat oktober de meest volatiele maand is geweest en december de minst volatiele. Belangrijk is dat we ook zien dat december in absolute zin de meer 'extreme' waarde is. Dit is belangrijk voor onze analyse in de volgende sectie.
#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()Mogelijk gemaakt door
Dat is onze data, nu door naar hypothesetoetsing…
Hypothesetoetsing: wat is de vraag?
Hypothesetoetsing is een van de meest fundamentele technieken in data science, maar ook een van de meest intimiderende en slecht begrepen. De basis van deze angst ligt in hoe het wordt onderwezen in Statistiek 101, waar we te horen krijgen:
"voer een t‑toets uit, is het een eenzijdige of tweezijdige toets? Kies een geschikte teststatistiek zoals Welch’ t‑toets, bereken vrijheidsgraden, bereken de t‑score, zoek de kritieke waarde op in een tabel, vergelijk de kritieke waarde met de t‑statistiek ……"
Begrijpelijk leidt dit tot verwarring over welke toets je moet uitvoeren en hoe je die uitvoert. Al deze klassieke statistische technieken voor het uitvoeren van een hypothesetoets zijn echter ontwikkeld in een tijd met weinig rekenkracht en zijn simpelweg gesloten‑vorm analytische oplossingen om een p‑waarde te berekenen, niets meer! Maar met de extra complicatie dat je de juiste formule moet kiezen voor de gegeven situatie, door hun restrictieve en soms ondoorzichtige aannames.
Maar verheug je!
Er is een betere manier. Simulatie.
Om te begrijpen hoe simulatie ons kan helpen, herinneren we ons wat een hypothesetoets is:
We willen toetsen "of het waargenomen effect in onze data echt is of dat het simpelweg toevallig kan gebeuren" en om deze toets uit te voeren doen we het volgende:
Kies een geschikte 'teststatistiek': dit is simpelweg een getal dat het waargenomen effect meet. In ons geval kiezen we de absolute afwijking in AMVR ten opzichte van het gemiddelde.
Construeer een nulhypothese: dit is een versie van de data waarin het waargenomen effect niet aanwezig is. In ons geval schudden we de labels van de data herhaaldelijk door elkaar (permutatie). De rechtvaardiging hiervoor staat hieronder uitgelegd.
Bereken een p‑waarde: dit is de kans om het waargenomen effect in de nuldata te zien, met andere woorden, door toeval. Dit doen we via herhaalde simulatie van de nuldata. In ons geval schudden we de 'date'‑labels van de data vaak en tellen we simpelweg hoe vaak onze teststatistiek voorkomt in meerdere simulaties.
Dat is hypothesetoetsing in 3 stappen! Ongeacht welk fenomeen we toetsen, de vraag is altijd dezelfde: "is het waargenomen effect echt, of is het toeval"
De echte kracht van simulatie is dat we expliciet moeten maken wat onze modelaannames zijn via code. Terwijl klassieke technieken een 'black‑box' kunnen zijn als het om hun aannames gaat.
Voorbeeld hieronder: De linkerplot toont de echte data en het waargenomen effect met een bepaalde waarschijnlijkheid (groen). De rechterplot is onze gesimuleerde nuldata met een registratie van wanneer het waargenomen effect toevallig werd gezien (rood). Dit is de basis van hypothesetoetsing: wat is de kans om het waargenomen effect te zien op onze nuldata.
Het meest kritieke deel van hypothesetoetsing is helder hebben welke vraag we proberen te beantwoorden. In ons geval vragen we:
"Zou de meest extreme waarde toevallig kunnen optreden?"
De meest extreme waarde definiëren we als de grootste absolute AMVR‑afwijking van het gemiddelde. Deze vraag vormt onze nulhypothese.
In onze data is de meest extreme waarde de decemberwaarde (1,23), niet de oktoberwaarde (1,08), omdat we kijken naar de meest significante absolute afwijking van het gemiddelde, niet simpelweg de hoogste volatiliteit.
1. Om de absolute afwijking te krijgen nemen we elk waarde‑minus‑gemiddelde en passen we de functie abs() toe.
2. Met de methode sort_values() kunnen we de resultaten van klein naar groot ordenen. Selecteer vervolgens de 2 grootste waarden, okt & dec.
#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)Mogelijk gemaakt door
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.080729166666667Mogelijk gemaakt door
Simulatie
Nu we weten welke vraag we stellen, moeten we ons 'nulmodel' construeren.
Er zijn een aantal opties:
Parametrische modellen. Als we een goed idee hadden van de dataverdeling of daar simpelweg aannames over deden, zouden we 'klassieke' hypothesetoetsingstechnieken kunnen gebruiken: t‑toets, X², one‑way ANOVA, enzovoort. Deze modellen kunnen beperkend zijn en enigszins een blackbox als de onderzoeker hun aannames niet volledig begrijpt.
Directe simulatie. We zouden aannames kunnen doen over het datagenererend proces en direct simuleren. Zo zouden we een ARMA‑tijdreeksmodel kunnen specificeren voor de financiële data die we modelleren en het bewust zo ontwerpen dat er geen seizoenspatroon is. Dit kan een redelijke keuze zijn voor ons probleem. Maar als we het datagenererend proces voor de S&P500 zouden kennen, waren we al rijk!
Simulatie via hersampling. Dit is de aanpak die we nemen. Door herhaaldelijk willekeurig te sampelen uit de bestaande dataset en de labels te schudden, kunnen we het waargenomen effect even waarschijnlijk maken onder alle labels in onze data (in ons geval zijn de labels de datums), en zo de gewenste nuldataset verkrijgen.
Sampling is een groot onderwerp, maar we richten ons op één specifieke techniek, permutatie of shufflen.
Om het gewenste nulmodel te krijgen, moeten we een dataset construeren waarin geen seizoenspatroon aanwezig is. Als de nul waar is, dat er geen seizoenspatroon in de data is en het waargenomen effect simpelweg toeval was, dan zijn de labels voor elke maand (jan, feb, enz.) betekenisloos en kunnen we de data herhaaldelijk schudden om op te bouwen wat klassieke statistiek de 'steekproefverdeling van de teststatistiek onder de nulhypothese' zou noemen. Dit heeft het gewenste effect dat het waargenomen fenomeen (de extreme decemberwaarde) even waarschijnlijk wordt voor alle maanden, precies wat ons nulmodel vereist.
Om te bewijzen hoe krachtig simulatie is met moderne rekenkracht, zal de code in dit voorbeeld daadwerkelijk de dagelijkse prijsdata permuteren, wat veel meer rekenwerk vereist, maar op een moderne CPU nog steeds in seconden voltooid is.
Let op: het schudden van ofwel de dagelijkse of de maandelijkse labels levert in ons geval de gewenste nuldataset op.
Een geweldige bron om over sampling te leren is van Julian Simon.
Let op: De manier waarop onze toets is geconstrueerd is equivalent aan een tweezijdige toets met 'klassieke' methoden, zoals Welch’ t‑toets of ANOVA, enz., omdat we geïnteresseerd zijn in de meest extreme waarde, boven of onder het gemiddelde.
Deze beslissingen zijn ontwerpkeuzes, en we hebben die vrijheid omdat het nulmodel precies dat is: een model! Dit betekent dat we de parameters kunnen specificeren zoals we willen; het belangrijkste is heel duidelijk te zijn over welke vraag we proberen te beantwoorden.
We willen dus Python gebruiken om veel data te simuleren om een nuldataset te creëren. Om dit te bereiken simuleren we 1000 sets van 12 AMVR, waarbij we elke keer de 'date'‑labels permuteren om de steekproefverdeling op te bouwen. De output van deze code is hieronder opgenomen in de p‑hacking‑sectie.
1. Eerst definiëren we een paar containers om onze resultaten op te slaan: we gebruiken een pd.DataFrame() en een simpele array [ ].
2. Vervolgens definiëren we een teller en zetten die op nul, en starten dan een for‑loop van 1000 iteraties om de gesimuleerde data te maken.
3. De eerste regel binnen de for‑loop neemt de oorspronkelijke dagelijkse rendementen (van het begin van de tutorial) en gebruikt de pandas‑methode sample(). Deze sampelt willekeurig uit de dagelijkse rendementen een gegeven aantal keren. In ons geval hebben we 8191 datapunten (252 handelsdagen per jaar over 32 jaar). We laten ook de index vallen met reset_index(), zodat we een nieuwe index kunnen toevoegen. Dit is nodig omdat de oorspronkelijke datumindex wordt mee geschud met de data, terwijl we alleen de data willen 'shufflen'.
4. De volgende regel voegt de nieuwe datumindex toe door de index van daily_ret_shuffle opnieuw toe te wijzen met een pd.bdate_range van gelijke lengte aan de oorspronkelijke data. We gebruiken bdate_range in plaats van date_range, omdat we werkdagen willen (5, ma–vr) in plaats van 7 dagen per week.
5. Nu we onze 'geschudde' data hebben gegenereerd, voeren we dezelfde datatransformaties uit als aan het begin om onze AMVR‑waarden te construeren. Dat is wat de volgende 3 regels code doen.
6. Nu we onze gesimuleerde AMVR‑waarden hebben, moeten we ze toevoegen aan onze dataframe die alle 1000 runs van de simulatie opslaat. pd.concat neemt onze huidige dataframe en voegt elke volgende dataframe eraan toe. We kiezen as 1, zodat we kolommen met data krijgen.
7. Maximonth is een variabele die alleen de hoogste waarde van elke simulatie opslaat (die we later in onze analyse gebruiken).
8. Omdat onze analyse absolute AMVR‑waarden vereist, nemen de volgende 3 regels code alle 1000 AMVR‑runs en vlakken die af tot één grote array. We nemen dan het gemiddelde en trekken de individuele waarden af van dit gemiddelde en nemen tenslotte de absolute waarde.
9. We doen hetzelfde voor de 'alleen hoogste' data. Let op: hier gebruikten we een lijst in plaats van een dataframe, wat betekent dat we een list comprehension kunnen gebruiken om de abs(AMVR) voor elk van onze hoogste waarden te berekenen.
We hebben nu alle data die we nodig hebben om onze analyse te voltooien: zowel de oorspronkelijke observaties als onze nieuwe gesimuleerde dataset.
#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]Mogelijk gemaakt door
p‑hacking
Hier wordt het interessant. We hebben een hypothese geconstrueerd om te toetsen, we hebben gesimuleerde data gegenereerd door de 'date'‑labels van de data te schudden, nu moeten we onze hypothesetoets uitvoeren om de kans te vinden om een resultaat zo significant als het decemberresultaat te observeren, gegeven dat de nulhypothese (geen seizoenspatroon) waar is.
Voordat we de toets uitvoeren, stellen we onze verwachtingen vast.
Wat is de kans om minstens één significant resultaat te zien gegeven een significantieniveau van 5%?
= 1‑p(niet significant)
= 1‑(1–0,05)¹²
= 0,46
er is dus 46% kans om minstens één maand met een significant resultaat te zien, gegeven dat onze nul waar is.
Vraag nu, voor elke individuele toets (waarbij we elke van de 12 maanden absolute AMVR vergelijken met het gemiddelde) hoeveel significante waarden verwachten we te zien in onze willekeurige, niet‑seizoensgebonden data?
12 x 0,05 = 0,6
Dus bij een significantieniveau van 0,05 mogen we een vals‑positieve rate van 0,6 verwachten. Met andere woorden, voor elke toets (met de nuldata) waarbij we alle 12 maanden AMVR vergelijken met het gemiddelde, zal 0,6 maand een significant resultaat laten zien. (we kunnen uiteraard niet minder dan 1 maand een resultaat laten zien, maar onder herhaalde toetsing zal de uitkomst naar dit getal neigen).
Gedurende dit werk hebben we het belang benadrukt van heel duidelijk zijn over de vraag die we proberen te beantwoorden. Het probleem met de verwachtingen die we zojuist hebben berekend, is dat we ervan uitgingen dat we toetsen op een significant resultaat tegen alle 12 maanden! Daarom is de kans om minstens één vals‑positief te zien zo hoog, 46%.
Dit is een voorbeeld van multiple comparison bias, waarbij we onze zoekruimte hebben vergroot en de kans om een significant resultaat te vinden is toegenomen. Dit is een probleem omdat we dit fenomeen zouden kunnen misbruiken om de parameters van ons model te cherry‑picken die ons de 'gewenste' p‑waarde geven.
Dit is de essentie van p‑hacking
Om het effect van p‑hacking te illustreren en hoe je multipliciteit vermindert, moeten we het subtiele maar belangrijke verschil begrijpen tussen de volgende 2 vragen:
"Wat is de kans dat december er zó extreem uitziet door toeval?"
"Wat is de kans dat om het even welke maand er zó extreem uitziet door toeval?"
De schoonheid van simulatie zit in de eenvoud. De volgende code is alles wat we nodig hebben om de p‑waarde te berekenen om de eerste vraag te beantwoorden. We tellen simpelweg hoeveel waarden in onze dataset met alle 12.000 AMVR‑afwijkingen (12 maanden x 1000 trials) groter zijn dan de geobserveerde decemberwaarde. We krijgen een p‑waarde van 4,4%, dicht bij onze enigszins willekeurige 5%‑grens, maar toch significant.
#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 )Mogelijk gemaakt door
Om de tweede vraag te beantwoorden en multipliciteit te vermijden, vergelijken we ons resultaat niet met de verdeling die is gemaakt met alle 12.000 AMVR‑afwijkingen, maar beschouwen we alleen de hoogste waarde uit elk van de 1000 trials met absolute AMVR. Dit geeft een p‑waarde van 23%, duidelijk niet significant!
#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 )Mogelijk gemaakt door
Nu we onze resultaten hebben, plotten we deze verdelingen om het effect van p‑hacking en de resultaten van onze analyse visueel te tonen:
1. Eerst gebruiken we np.quantile() om ons 5%‑significantieniveau te vinden; we plotten dit met axvline in de onderste grafieken.
2. Vervolgens definiëren we onze 4 subplots. Hier retourneert plt.subplots 'fig' = onze figuur en 'ax1,ax2,ax3,ax4' = variabelen die elk van de 4 subplots vertegenwoordigen.
3. #plot 1 toont de eerste kolom. We definiëren eenvoudig een histogram van onze data 'abs_all_months' en kiezen type='bar'. Voor de onderste grafieken zetten we cumulative='True', dit geeft ons een cdf in plaats van een pdf. Bins bepaalt hoeveel balken we hebben om de waarden te vullen; 30 is een redelijk aantal voor onze data. We formatteren de grafieken vervolgens met technieken die we al hebben geleerd, zoals axvline om de significantie‑ en observatiewaarden te plotten.
abs_all_months_95 = np.quantile(abs_all_months,.95)
abs_highest_95 = np.quantile(abs_highest,.95)
fig,((ax1,ax2),(ax3,ax4))= plt.subplots(2,2,sharex='col',figsize=(20,20))#plot 1
ax1.hist(abs_all_months,histtype='bar',color='#42a5f5')
ax1.set_title('AMVR all months',fontsize=30)
ax1.set_ylabel('Frequency',fontsize=20)
ax3.hist(abs_all_months,density=1,histtype='bar',cumulative=True,bins=30,color='#42a5f5')
ax3.set_ylabel('Cumulative probability',fontsize=20)
ax1.axvline(Dec_value,color='b',label='Dec Result',lw=10)
ax3.axvline(Dec_value,color='b',lw=10)
ax3.axvline(abs_all_months_95,color='r',ls='--',label='5% Sig level',lw=10)#plot2
ax2.hist(abs_highest,histtype='bar',color='g')
ax2.set_title('AMVR highest only',fontsize=30)
ax2.axvline(Dec_value,color='b',lw=10)
ax4.hist(abs_highest,density=1,histtype='bar',cumulative=True,bins=30,color='g')
ax4.axvline(Dec_value,color='b',lw=10)
ax4.axvline(abs_highest_95,color='r',ls='--',lw=10)
ax1.legend(fontsize=15)
ax3.legend(fontsize=15)Mogelijk gemaakt door
De linker kolom bevat de data die vraag 1 beantwoordt en de rechter kolom vraag 2. De bovenste rij zijn de kansverdelingen, en de onderste rij is de CDF. De rode stippellijn is het 5%‑significantieniveau dat we willekeurig kozen. De blauwe lijn is de oorspronkelijke extreme december‑AMVR‑waarde van 1,23.
De plot aan de linkerkant laat zien dat de oorspronkelijke decemberwaarde significant is op 5%‑niveau, maar nipt! Wanneer we echter rekening houden met multiple comparison bias, verschuift in de rechterplot de drempel voor significantie omhoog van rond 1,2 (abs AMVR) naar ongeveer 1,6 (zie de rode lijn).
Door te corrigeren voor multiple comparison bias is onze decemberwaarde van 1,23 niet langer significant!
Door rekening te houden met de specifieke vraag die we proberen te beantwoorden en multiple comparison bias te vermijden, hebben we p‑hacking van ons model voorkomen en voorkomen dat we een significant resultaat tonen wanneer dat er niet is.
Wil je p‑hacking verder verkennen en zien hoe het kan worden misbruikt om een bepaald verhaal over onze data te vertellen? Bekijk dan deze geweldige interactieve app van FiveThirtyEight
Conclusies
We hebben geleerd dat hypothesetoetsing niet het enge monster is dat we dachten. Volg simpelweg de 3 stappen hierboven om je model te construeren voor elk soort data of teststatistiek.
We hebben laten zien dat de juiste vraag stellen cruciaal is voor wetenschappelijke analyse. Een kleine wijziging in de formulering kan leiden tot een heel ander model met heel andere resultaten.
Hopelijk is de kracht van de wat meer intermediate Python‑functionaliteit en z’n vermogen om je statistische toetsen te laten uitvoeren nu duidelijk! Met slechts een paar regels code hebben we een reëel fenomeen geconstrueerd en getoetst en konden we actiegerichte conclusies trekken.
We bespraken het belang van het herkennen en corrigeren voor multiple comparison bias en het vermijden van de valkuilen van p‑hacking, en lieten zien hoe een ogenschijnlijk significant resultaat niet‑significant kan worden.
Met steeds meer 'big data', samen met academische druk om een paper met 'nieuwe' bevindingen te produceren of politieke druk om een resultaat als 'significant' te tonen, neemt de verleiding tot p‑hacking toe. Door te leren herkennen wanneer we ons eraan bezondigen en daar passend voor te corrigeren, kunnen we betere onderzoekers worden en uiteindelijk nauwkeuriger en daarmee actiegerichte wetenschappelijke resultaten produceren!
Aantekeningen van de auteur: Onze resultaten verschillen enigszins van het oorspronkelijke Winton‑onderzoek, deels doordat we een enigszins andere dataset hebben (32 jaar vs. 87 jaar) en zij oktober als interessante maand hebben terwijl wij december hebben. Ook gebruikten zij een niet‑openbaar gemaakte methode voor hun 'gesimuleerde data', terwijl wij onze methodologie expliciet hebben gemaakt via code. We hebben gedurende dit werk bepaalde modelaannames gedaan; ook die zijn expliciet gemaakt en in de code te zien. Deze ontwerp‑ en modelkeuzes maken deel uit van het wetenschappelijke proces; zolang ze expliciet worden gemaakt, heeft de analyse waarde.
Bouw je Python-skills uit voor je carrière in finance. Leer werken met lijsten, arrays en datavisualisaties om financiële analyses onder de knie te krijgen.
%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)
Adel Nehme · 15 min
"In de afgelopen 32 jaar was oktober gemiddeld de meest volatiele maand voor de S&P500 en december de minst volatiele".
In deze tutorial gebruiken we Python om een volledige analyse en toetsing van dit fenomeen uit te voeren en te bepalen of het statistisch significant is of niet.
We gebruiken tijdreeksfunctionaliteit van pandas, zoals Resample, om ruwe aandelenkoersdata te bewerken tot een formaat dat de data tot leven brengt.
Vervolgens ontdekken we hoe je hypothesetoetsing uitvoert met simulatie in Python in plaats van omslachtige formules.
Tot slot bespreken we een groot probleem in statistische analyse, Multiple Comparison Bias, leren we hoe je daarmee omgaat en laten we via matplotlib visueel het effect van 'p-hacking' zien.
Ons doel:
Onze data:
We gebruiken dagelijkse S&P500-data voor deze analyse; specifiek gebruiken we de ruwe dagelijkse slotkoersen van 1986 tot 2018 (die verrassend moeilijk te vinden zijn, dus ik heb ze openbaar beschikbaar gemaakt).
De inspiratie voor dit artikel kwam van Winton, wat we hier zullen reproduceren, zij het met 32 jaar data tegenover hun 87 jaar.
Aan de slag met Pandas
Om de vraag te beantwoorden of de extreme volatiliteit in bepaalde maanden echt significant is en dus waarschijnlijk zal aanhouden, moeten we onze 32 jaar aan prijsdata omzetten naar een formaat dat het fenomeen dat we onderzoeken zichtbaar maakt.
Ons voorkeursformaat is de gemiddelde maandelijkse volatiliteitsrangschikking (AMVR).
De volgende code laat zien hoe we onze ruwe prijsdata in dit formaat krijgen. Aan de slag!
Eerst de standaard imports. (matplotlib.patches geeft ons controle over de styling van individuele balken in een histogram)
Een handige truc om grafieken op volledige breedte weer te geven in Jupyter-notebooks:
Deze is optioneel - probeer het eens en kijk wat er gebeurt!
Importeer onze data met de methode read_csv, die een pad krijgt, in ons geval een URL naar de data. We geven ook aan dat 'date' als index moet worden gebruikt en dat de datums automatisch zo goed mogelijk uit de tekst moeten worden geparseerd.
Dit levert de onbewerkte, niet-aangepaste slotkoersen van de S&P500 (SPX) op. Nu moeten we onze ruwe prijzen omzetten in dagelijkse %‑rendementen. Daarvoor hebben we twee opties:
Voor ons doel gebruiken we methode 2, omdat pandas de berekening direct aankan op een dataset van deze omvang (meer dan 8000 waarden). We schonen de data ook op door de eerste waarde te verwijderen, die 'NaN' wordt omdat er geen prijsverandering is ten opzichte van de dag ervoor. pct_change() neemt een optionele parameter 'periods' die de periodeshift verandert. We laten die op de standaardwaarde één staan.
De volgende stap is om deze dagelijkse %‑veranderingen te transformeren naar 'maandelijkse geannualiseerde volatiliteit'. De eerste regel code hieronder laat zien dat we deze transformatie in één regel kunnen uitvoeren. Dit is de kracht van pandas! Maar laten we elke stap op volgorde uitpakken.
1. Om mnthly_annu te krijgen, gebruiken we eerst de methode 'resample' op onze dagelijkse rendementen. Resample laat ons de frequentie van de dataperiodes veranderen. Het neemt een 'frequency offset string' of, simpel gezegd, een letter die overeenkomt met de gewenste nieuwe frequentie. Deze omvatten:
Aangezien we willen resamplen van dagelijkse rendementen naar maandelijks, geven we 'M' door als parameter.
2. Ten tweede moeten we beslissen hoe we tot deze nieuwe maandwaarde willen komen; we zouden de dagelijkse waarden kunnen sommeren, vermenigvuldigen, enzovoort. Voor onze analyse willen we een maat voor volatiliteit en standaarddeviatie werkt goed, dus voegen we std() toe aan onze geresamplede data om maandelijkse volatiliteit te krijgen.
3. De laatste stap is om dit cijfer te annualiseren. Dat doen we eenvoudig door te vermenigvuldigen met de wortel uit 12, waarbij 12 het aantal periodes (maanden) per jaar is. Dit levert de geannualiseerde maandvolatiliteiten op die we nodig hebben.
Best practice in data‑analyse is om de data te visualiseren terwijl we door de analyse heen gaan. Laten we dus naar onze geannualiseerde volatiliteitsdata kijken.
De onderstaande plot toont duidelijk grote marktgebeurtenissen, zoals Black Monday en de financiële crisis van 2008. De matplotlib‑methode 'axvspan' laat ons verticale lijnen toevoegen (axhspan is de overeenkomstige methode voor horizontale lijnen). Hij neemt 'xmin' en 'xmax' als parameters, die de breedte van de gecreëerde rechthoek specificeren; omdat onze x‑as de datetime‑index is, kunnen we eenvoudig de jaren doorgeven die we willen markeren. De alpha‑parameter past de transparantie van de blokkleur aan zodat we de grafiek eronder nog kunnen zien.
Het attribuut mpatches laat ons een aangepaste legenda maken. Dat doen we door eerst de variabele 'labs' te definiëren met kleur, alpha en de gewenste tekst, en die vervolgens door te geven aan de 'handles'‑parameter binnen plt.legend om de legenda te renderen.
We hebben nu een van de krachtige methodes in pandas gezien: resample. Laten we nu een andere gebruiken, 'groupby'. We moeten van onze maandelijkse geannualiseerde volatiliteitswaarden naar onze gewenste AMVR‑maatstaf. Ook dit kunnen we in slechts een paar regels code bereiken.
1. Eerst passen we de 'groupby'‑methode toe op onze mnthly_annu‑waarden. Groupby vereist een parameter die specificeert hoe de data gegroepeerd moet worden. Dit kan een functie zijn of, zoals in ons geval, een Series. We geven 'mnthly_annu.index.year' door; dit is simpelweg het jaarattribuut op de datetime‑index van de mnthly_annu‑data. Dit groepeert onze maandvolatiliteitswaarden voor elk van de 32 jaar in onze dataset.
2. Vervolgens passen we de methode 'rank' toe, die de data in oplopende volgorde rangschikt.
3. Tot slot doen we hetzelfde nog eens en middelen we over alle jaren per maand. Dit levert onze definitieve AMVR‑waarden op!
Dit levert onze definitieve Gemiddelde Maandelijkse Volatiliteitsrangschikking op. Numeriek zien we dat maand 10 (oktober) het hoogst is en 12 (december) het laagst.
De juiste visualisatie kiezen voor onze data is belangrijk, omdat we het verhaal van de data willen vertellen met de meest impactvolle grafieken en plots. Voor onze analyse willen we duidelijk de hoogste en de laagste volatiliteit tonen. Daarvoor gebruiken we de staafdiagram van matplotlib en annoteren en kleuren we dienovereenkomstig voor het meeste effect.
1. Door b_plot (b_plot[9]) te indexeren met de gewenste balk kunnen we de kleur instellen om de hoogste en laagste waarden te markeren.
2. Om de AMVR‑waarden aan hun respectieve balk toe te voegen, lopen we door elke AMVR‑waarde ('final'‑variabele) en ronden we het resultaat af op 2 decimalen. 'i' en 'v' staan voor de 'index' en 'waarde' van de 'final'‑variabele, dus 'i' loopt van 1–12 en 'value' geeft de AMVR‑waarde afgerond op 2 decimalen. We gebruiken deze variabelen binnen plt.text(). De eerste parameter van plt.txt() is de x‑aslocatie voor het label, dus we geven 'i' door voor elke balk en verschuiven met 0,8 om te centreren. De tweede parameter is de stringversie van de 'final'‑waarde.
3. Om de gemiddelde waarde te tonen gebruiken we axhline en gebruiken we de 'label'‑waarde om rechtsboven een legenda toe te voegen. Geef simpelweg de y‑aswaarde door, oftewel het gemiddelde, en voeg vervolgens styling toe.
Dit levert de beste visual tot nu toe van het fenomeen dat we onderzoeken. We zien duidelijk dat oktober de meest volatiele maand is geweest en december de minst volatiele. Belangrijk is dat we ook zien dat december in absolute zin de meer 'extreme' waarde is. Dit is belangrijk voor onze analyse in de volgende sectie.
Dat is onze data, nu door naar hypothesetoetsing…
Hypothesetoetsing: wat is de vraag?
Hypothesetoetsing is een van de meest fundamentele technieken in data science, maar ook een van de meest intimiderende en slecht begrepen. De basis van deze angst ligt in hoe het wordt onderwezen in Statistiek 101, waar we te horen krijgen:
"voer een t‑toets uit, is het een eenzijdige of tweezijdige toets? Kies een geschikte teststatistiek zoals Welch’ t‑toets, bereken vrijheidsgraden, bereken de t‑score, zoek de kritieke waarde op in een tabel, vergelijk de kritieke waarde met de t‑statistiek ……"
Begrijpelijk leidt dit tot verwarring over welke toets je moet uitvoeren en hoe je die uitvoert. Al deze klassieke statistische technieken voor het uitvoeren van een hypothesetoets zijn echter ontwikkeld in een tijd met weinig rekenkracht en zijn simpelweg gesloten‑vorm analytische oplossingen om een p‑waarde te berekenen, niets meer! Maar met de extra complicatie dat je de juiste formule moet kiezen voor de gegeven situatie, door hun restrictieve en soms ondoorzichtige aannames.
Maar verheug je!
Er is een betere manier. Simulatie.
Om te begrijpen hoe simulatie ons kan helpen, herinneren we ons wat een hypothesetoets is:
We willen toetsen "of het waargenomen effect in onze data echt is of dat het simpelweg toevallig kan gebeuren" en om deze toets uit te voeren doen we het volgende:
Er is maar één toets! Deze geweldige blog van Allen Downey geeft meer details over hypothesetoetsing
De echte kracht van simulatie is dat we expliciet moeten maken wat onze modelaannames zijn via code. Terwijl klassieke technieken een 'black‑box' kunnen zijn als het om hun aannames gaat.
Voorbeeld hieronder: De linkerplot toont de echte data en het waargenomen effect met een bepaalde waarschijnlijkheid (groen). De rechterplot is onze gesimuleerde nuldata met een registratie van wanneer het waargenomen effect toevallig werd gezien (rood). Dit is de basis van hypothesetoetsing: wat is de kans om het waargenomen effect te zien op onze nuldata.
Het meest kritieke deel van hypothesetoetsing is helder hebben welke vraag we proberen te beantwoorden. In ons geval vragen we:
"Zou de meest extreme waarde toevallig kunnen optreden?"
De meest extreme waarde definiëren we als de grootste absolute AMVR‑afwijking van het gemiddelde. Deze vraag vormt onze nulhypothese.
In onze data is de meest extreme waarde de decemberwaarde (1,23), niet de oktoberwaarde (1,08), omdat we kijken naar de meest significante absolute afwijking van het gemiddelde, niet simpelweg de hoogste volatiliteit.
1. Om de absolute afwijking te krijgen nemen we elk waarde‑minus‑gemiddelde en passen we de functie abs() toe.
2. Met de methode sort_values() kunnen we de resultaten van klein naar groot ordenen. Selecteer vervolgens de 2 grootste waarden, okt & dec.
Simulatie
Nu we weten welke vraag we stellen, moeten we ons 'nulmodel' construeren.
Er zijn een aantal opties:
Sampling is een groot onderwerp, maar we richten ons op één specifieke techniek, permutatie of shufflen.
Om het gewenste nulmodel te krijgen, moeten we een dataset construeren waarin geen seizoenspatroon aanwezig is. Als de nul waar is, dat er geen seizoenspatroon in de data is en het waargenomen effect simpelweg toeval was, dan zijn de labels voor elke maand (jan, feb, enz.) betekenisloos en kunnen we de data herhaaldelijk schudden om op te bouwen wat klassieke statistiek de 'steekproefverdeling van de teststatistiek onder de nulhypothese' zou noemen. Dit heeft het gewenste effect dat het waargenomen fenomeen (de extreme decemberwaarde) even waarschijnlijk wordt voor alle maanden, precies wat ons nulmodel vereist.
Om te bewijzen hoe krachtig simulatie is met moderne rekenkracht, zal de code in dit voorbeeld daadwerkelijk de dagelijkse prijsdata permuteren, wat veel meer rekenwerk vereist, maar op een moderne CPU nog steeds in seconden voltooid is.
Let op: het schudden van ofwel de dagelijkse of de maandelijkse labels levert in ons geval de gewenste nuldataset op.
Een geweldige bron om over sampling te leren is van Julian Simon.
Let op: De manier waarop onze toets is geconstrueerd is equivalent aan een tweezijdige toets met 'klassieke' methoden, zoals Welch’ t‑toets of ANOVA, enz., omdat we geïnteresseerd zijn in de meest extreme waarde, boven of onder het gemiddelde.
Deze beslissingen zijn ontwerpkeuzes, en we hebben die vrijheid omdat het nulmodel precies dat is: een model! Dit betekent dat we de parameters kunnen specificeren zoals we willen; het belangrijkste is heel duidelijk te zijn over welke vraag we proberen te beantwoorden.
We willen dus Python gebruiken om veel data te simuleren om een nuldataset te creëren. Om dit te bereiken simuleren we 1000 sets van 12 AMVR, waarbij we elke keer de 'date'‑labels permuteren om de steekproefverdeling op te bouwen. De output van deze code is hieronder opgenomen in de p‑hacking‑sectie.
1. Eerst definiëren we een paar containers om onze resultaten op te slaan: we gebruiken een pd.DataFrame() en een simpele array [ ].
2. Vervolgens definiëren we een teller en zetten die op nul, en starten dan een for‑loop van 1000 iteraties om de gesimuleerde data te maken.
3. De eerste regel binnen de for‑loop neemt de oorspronkelijke dagelijkse rendementen (van het begin van de tutorial) en gebruikt de pandas‑methode sample(). Deze sampelt willekeurig uit de dagelijkse rendementen een gegeven aantal keren. In ons geval hebben we 8191 datapunten (252 handelsdagen per jaar over 32 jaar). We laten ook de index vallen met reset_index(), zodat we een nieuwe index kunnen toevoegen. Dit is nodig omdat de oorspronkelijke datumindex wordt mee geschud met de data, terwijl we alleen de data willen 'shufflen'.
4. De volgende regel voegt de nieuwe datumindex toe door de index van daily_ret_shuffle opnieuw toe te wijzen met een pd.bdate_range van gelijke lengte aan de oorspronkelijke data. We gebruiken bdate_range in plaats van date_range, omdat we werkdagen willen (5, ma–vr) in plaats van 7 dagen per week.
5. Nu we onze 'geschudde' data hebben gegenereerd, voeren we dezelfde datatransformaties uit als aan het begin om onze AMVR‑waarden te construeren. Dat is wat de volgende 3 regels code doen.
6. Nu we onze gesimuleerde AMVR‑waarden hebben, moeten we ze toevoegen aan onze dataframe die alle 1000 runs van de simulatie opslaat. pd.concat neemt onze huidige dataframe en voegt elke volgende dataframe eraan toe. We kiezen as 1, zodat we kolommen met data krijgen.
7. Maximonth is een variabele die alleen de hoogste waarde van elke simulatie opslaat (die we later in onze analyse gebruiken).
8. Omdat onze analyse absolute AMVR‑waarden vereist, nemen de volgende 3 regels code alle 1000 AMVR‑runs en vlakken die af tot één grote array. We nemen dan het gemiddelde en trekken de individuele waarden af van dit gemiddelde en nemen tenslotte de absolute waarde.
9. We doen hetzelfde voor de 'alleen hoogste' data. Let op: hier gebruikten we een lijst in plaats van een dataframe, wat betekent dat we een list comprehension kunnen gebruiken om de abs(AMVR) voor elk van onze hoogste waarden te berekenen.
We hebben nu alle data die we nodig hebben om onze analyse te voltooien: zowel de oorspronkelijke observaties als onze nieuwe gesimuleerde dataset.
p‑hacking
Hier wordt het interessant. We hebben een hypothese geconstrueerd om te toetsen, we hebben gesimuleerde data gegenereerd door de 'date'‑labels van de data te schudden, nu moeten we onze hypothesetoets uitvoeren om de kans te vinden om een resultaat zo significant als het decemberresultaat te observeren, gegeven dat de nulhypothese (geen seizoenspatroon) waar is.
Voordat we de toets uitvoeren, stellen we onze verwachtingen vast.
Wat is de kans om minstens één significant resultaat te zien gegeven een significantieniveau van 5%?
= 1‑p(niet significant)
= 1‑(1–0,05)¹²
= 0,46
er is dus 46% kans om minstens één maand met een significant resultaat te zien, gegeven dat onze nul waar is.
Vraag nu, voor elke individuele toets (waarbij we elke van de 12 maanden absolute AMVR vergelijken met het gemiddelde) hoeveel significante waarden verwachten we te zien in onze willekeurige, niet‑seizoensgebonden data?
12 x 0,05 = 0,6
Dus bij een significantieniveau van 0,05 mogen we een vals‑positieve rate van 0,6 verwachten. Met andere woorden, voor elke toets (met de nuldata) waarbij we alle 12 maanden AMVR vergelijken met het gemiddelde, zal 0,6 maand een significant resultaat laten zien. (we kunnen uiteraard niet minder dan 1 maand een resultaat laten zien, maar onder herhaalde toetsing zal de uitkomst naar dit getal neigen).
Gedurende dit werk hebben we het belang benadrukt van heel duidelijk zijn over de vraag die we proberen te beantwoorden. Het probleem met de verwachtingen die we zojuist hebben berekend, is dat we ervan uitgingen dat we toetsen op een significant resultaat tegen alle 12 maanden! Daarom is de kans om minstens één vals‑positief te zien zo hoog, 46%.
Dit is een voorbeeld van multiple comparison bias, waarbij we onze zoekruimte hebben vergroot en de kans om een significant resultaat te vinden is toegenomen. Dit is een probleem omdat we dit fenomeen zouden kunnen misbruiken om de parameters van ons model te cherry‑picken die ons de 'gewenste' p‑waarde geven.
Om het effect van p‑hacking te illustreren en hoe je multipliciteit vermindert, moeten we het subtiele maar belangrijke verschil begrijpen tussen de volgende 2 vragen:
De schoonheid van simulatie zit in de eenvoud. De volgende code is alles wat we nodig hebben om de p‑waarde te berekenen om de eerste vraag te beantwoorden. We tellen simpelweg hoeveel waarden in onze dataset met alle 12.000 AMVR‑afwijkingen (12 maanden x 1000 trials) groter zijn dan de geobserveerde decemberwaarde. We krijgen een p‑waarde van 4,4%, dicht bij onze enigszins willekeurige 5%‑grens, maar toch significant.
Om de tweede vraag te beantwoorden en multipliciteit te vermijden, vergelijken we ons resultaat niet met de verdeling die is gemaakt met alle 12.000 AMVR‑afwijkingen, maar beschouwen we alleen de hoogste waarde uit elk van de 1000 trials met absolute AMVR. Dit geeft een p‑waarde van 23%, duidelijk niet significant!
Nu we onze resultaten hebben, plotten we deze verdelingen om het effect van p‑hacking en de resultaten van onze analyse visueel te tonen:
1. Eerst gebruiken we np.quantile() om ons 5%‑significantieniveau te vinden; we plotten dit met axvline in de onderste grafieken.
2. Vervolgens definiëren we onze 4 subplots. Hier retourneert plt.subplots 'fig' = onze figuur en 'ax1,ax2,ax3,ax4' = variabelen die elk van de 4 subplots vertegenwoordigen.
3. #plot 1 toont de eerste kolom. We definiëren eenvoudig een histogram van onze data 'abs_all_months' en kiezen type='bar'. Voor de onderste grafieken zetten we cumulative='True', dit geeft ons een cdf in plaats van een pdf. Bins bepaalt hoeveel balken we hebben om de waarden te vullen; 30 is een redelijk aantal voor onze data. We formatteren de grafieken vervolgens met technieken die we al hebben geleerd, zoals axvline om de significantie‑ en observatiewaarden te plotten.
De linker kolom bevat de data die vraag 1 beantwoordt en de rechter kolom vraag 2. De bovenste rij zijn de kansverdelingen, en de onderste rij is de CDF. De rode stippellijn is het 5%‑significantieniveau dat we willekeurig kozen. De blauwe lijn is de oorspronkelijke extreme december‑AMVR‑waarde van 1,23.
De plot aan de linkerkant laat zien dat de oorspronkelijke decemberwaarde significant is op 5%‑niveau, maar nipt! Wanneer we echter rekening houden met multiple comparison bias, verschuift in de rechterplot de drempel voor significantie omhoog van rond 1,2 (abs AMVR) naar ongeveer 1,6 (zie de rode lijn).
Door rekening te houden met de specifieke vraag die we proberen te beantwoorden en multiple comparison bias te vermijden, hebben we p‑hacking van ons model voorkomen en voorkomen dat we een significant resultaat tonen wanneer dat er niet is.
Wil je p‑hacking verder verkennen en zien hoe het kan worden misbruikt om een bepaald verhaal over onze data te vertellen? Bekijk dan deze geweldige interactieve app van FiveThirtyEight
Conclusies
Aantekeningen van de auteur: Onze resultaten verschillen enigszins van het oorspronkelijke Winton‑onderzoek, deels doordat we een enigszins andere dataset hebben (32 jaar vs. 87 jaar) en zij oktober als interessante maand hebben terwijl wij december hebben. Ook gebruikten zij een niet‑openbaar gemaakte methode voor hun 'gesimuleerde data', terwijl wij onze methodologie expliciet hebben gemaakt via code. We hebben gedurende dit werk bepaalde modelaannames gedaan; ook die zijn expliciet gemaakt en in de code te zien. Deze ontwerp‑ en modelkeuzes maken deel uit van het wetenschappelijke proces; zolang ze expliciet worden gemaakt, heeft de analyse waarde.
Als je meer wilt leren over finance in Python, volg dan DataCamp’s Intro to Python for Finance‑cursus en bekijk de Python For Finance Tutorial: Algorithmic Trading.
Volg me op twitter.com/pdquant voor meer!