Saham, Uji Signifikansi & p-Hacking

Oktober secara historis merupakan bulan paling volatil untuk saham, tetapi apakah ini sinyal yang persisten atau sekadar derau dalam data?

"Selama 32 tahun terakhir, Oktober rata-rata merupakan bulan paling volatil untuk S&P500 dan Desember yang paling rendah volatilitasnya".

Dalam tutorial ini, kita akan menggunakan Python untuk menelusuri analisis dan pengujian menyeluruh terhadap fenomena ini guna memastikan apakah secara statistik signifikan atau tidak.

Kita akan memanfaatkan fungsionalitas deret waktu pandas, seperti Resample, untuk mengolah data harga saham mentah ke format yang membuat data lebih hidup.

Kita kemudian akan mempelajari cara melakukan pengujian hipotesis menggunakan simulasi di Python alih-alih rumus yang merepotkan.

Terakhir, kita akan membahas masalah besar dalam analisis statistik, Multiple Comparison Bias, mempelajari cara menanganinya dan menampilkan secara visual melalui matplotlib, dampak dari 'p-Hacking'.

Tujuan kami:

• Menunjukkan cara menggunakan Pandas untuk menganalisis Deret Waktu
• Memahami cara membangun uji hipotesis
• Menggunakan simulasi dengan Python untuk melakukan uji hipotesis
• Menunjukkan pentingnya memperhitungkan bias perbandingan ganda (multiple comparison)

Data kami:

Kita akan menggunakan data harian S&P500 untuk analisis ini, khususnya, kita akan menggunakan harga penutupan harian mentah dari 1986 hingga 2018 (yang ternyata sulit dicari, jadi saya telah membuatnya tersedia untuk publik).

Inspirasi untuk tulisan ini datang dari Winton, yang akan kita replikasi di sini, meski dengan 32 tahun data vs. 87 tahun milik mereka.

Mengolah data dengan Pandas

Untuk menjawab apakah volatilitas ekstrem pada bulan-bulan tertentu benar-benar signifikan dan karenanya cenderung berlanjut, kita perlu mentransformasikan 32 tahun data harga menjadi format yang menampilkan fenomena yang sedang kita selidiki.

Format pilihan kita adalah peringkat volatilitas bulanan rata-rata (Average Monthly Volatility Rankings/AMVR).

Kode berikut menunjukkan cara kita mengubah data harga mentah ke format ini. Mari mulai!

Pertama impor standar. (matplotlib.patches memberi kita kontrol atas gaya setiap batang dalam histogram)

#standard imports
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
%matplotlib inline
Didukung oleh

Trik praktis agar bagan tampil selebar layar di notebook Jupyter:

#resize charts to fit screen if using Jupyter Notebook
plt.rcParams['figure.figsize']=[15,5]
Didukung oleh

Yang ini opsional — coba saja dan lihat apa yang terjadi!

#plt.xkcd()
Didukung oleh

Impor data kita menggunakan metode read_csv yang menerima path, dalam kasus kita berupa URL ke data. Kita juga mengatur agar menggunakan 'date' sebagai indeks dan secara otomatis mengurai tanggal dari teks yang diberikan, sebisa mungkin.

#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()
Didukung oleh

Ini memberi kita harga penutupan mentah yang belum disesuaikan dari S&P500 (SPX). Sekarang kita perlu mengonversi harga mentah menjadi imbal hasil harian dalam %. Untuk ini, kita punya dua opsi:

1. Kita bisa mengambil log natural dari harga. Ini akan memberi pendekatan terhadap imbal hasil harian yang sebenarnya.
2. Kita bisa menggunakan metode pandas 'pct_change()' untuk menghitung perubahan persentase harian secara langsung.

Untuk tujuan kita, kita akan menggunakan metode 2, karena pandas dapat menangani komputasi secara instan pada dataset sebesar ini (lebih dari 8000 nilai). Kita juga akan membersihkan data dengan membuang nilai pertama yang menjadi 'NaN' karena tidak ada perubahan harga dari hari sebelumnya. pct_change() memiliki parameter opsional 'periods' yang mengubah pergeseran periode. Kita biarkan pada default satu.

#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()
Didukung oleh

date
1986-01-03    0.006155
1986-01-06   -0.001091
1986-01-07    0.014954
1986-01-08   -0.027268
1986-01-09   -0.008944
Name: close, dtype: float64
Didukung oleh

Langkah berikutnya adalah mengambil nilai perubahan % harian ini dan mengubahnya menjadi 'volatilitas bulanan yang dianualkan'. Baris kode pertama di bawah menunjukkan bahwa kita bisa melakukan transformasi ini hanya dalam satu baris kode. Inilah kekuatan Pandas! Namun mari kita uraikan setiap langkahnya secara berurutan.

 1. Untuk mendapatkan mnthly_annu, pertama kita gunakan metode 'resample' pada imbal hasil harian. Resample memungkinkan kita mengubah frekuensi periode data. Ia menerima 'frequency offset string' atau sederhananya, huruf yang sesuai dengan frekuensi baru yang diinginkan. Di antaranya:

• B frekuensi hari kerja
• C frekuensi hari kerja kustom
• D frekuensi hari kalender
• W frekuensi mingguan
• M frekuensi akhir bulan

Karena kita ingin mengubah sampel dari imbal hasil harian ke bulanan, kita berikan 'M' sebagai parameter.

 2. Hal kedua yang perlu kita lakukan adalah memutuskan bagaimana kita ingin mendapatkan angka bulanan ini, misalnya, kita bisa menjumlahkan nilai harian, mengalikan, dll. Untuk analisis kita, kita ingin ukuran volatilitas dan deviasi standar bekerja dengan baik, jadi kita tambahkan std() ke data yang di-resample untuk mendapatkan volatilitas bulanan.

 3. Langkah terakhir adalah menganualkan angka ini. Kita lakukan dengan mengalikan akar kuadrat dari 12, dengan 12 adalah jumlah periode (bulan) dalam setahun. Ini memberi kita nilai volatilitas bulanan yang dianualkan yang kita butuhkan.

Praktik terbaik dalam analisis data adalah memvisualisasikan data saat kita melangkah. Jadi mari lihat data volatilitas yang dianualkan kita.

Plot di bawah menunjukkan peristiwa pasar besar dengan jelas, seperti Black Monday dan Krisis Keuangan 2008. Metode matplotlib 'axvspan' memungkinkan kita menambahkan garis vertikal (axhspan adalah metode terkait untuk garis horizontal), ia menerima 'xmin' dan 'xmax' sebagai parameter, yang menentukan lebar persegi panjang yang dibuat; karena sumbu x kita adalah indeks datetime, kita cukup meneruskan tahun yang ingin kita sorot. Parameter alpha menyesuaikan transparansi warna blok agar grafik di bawahnya tetap terlihat.

Atribut mpatches memungkinkan kita membuat kotak legenda kustom. Kita lakukan dengan terlebih dahulu mendefinisikan variabel 'labs' dengan warna, alpha, dan teks yang diinginkan, lalu meneruskannya ke parameter 'handles' dalam plt.legend untuk menampilkan 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])
Didukung oleh

date
1986-01-31    0.033317
1986-02-28    0.023585
1986-03-31    0.027961
1986-04-30    0.037426
1986-05-31    0.027412
Freq: M, Name: close, dtype: float64

<matplotlib.legend.Legend at 0x280b1ee6908>
Didukung oleh

Jadi kita telah melihat salah satu metode kuat di pandas, resample. Sekarang, mari gunakan metode lain, 'groupby'. Kita perlu beralih dari nilai volatilitas bulanan yang dianualkan ke metrik AMVR yang kita inginkan. Lagi-lagi, kita bisa mencapainya hanya dengan beberapa baris kode.

 1. Pertama, kita terapkan metode 'groupby' pada nilai mnthly_annu. Groupby memerlukan parameter yang menentukan bagaimana mengelompokkan data. Ini bisa berupa fungsi atau seperti kasus kita, sebuah Series. Kita teruskan 'mnthly_annu.index.year', ini hanyalah atribut tahun pada indeks datetime dari data mnthly_annu. Ini akan mengelompokkan nilai volatilitas bulanan untuk masing-masing dari 32 tahun dalam dataset kita.

 2. Selanjutnya, kita terapkan metode 'rank', yang mengurutkan data secara menaik.

 3. Terakhir, kita lakukan hal yang sama lagi dan rata-ratakan di seluruh tahun untuk setiap bulan. Ini memberi kita nilai AMVR akhir!

#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()
Didukung oleh

count    12.000000
mean      6.450521
std       0.627458
min       5.218750
25%       6.031013
50%       6.491004
75%       6.704545
max       7.531250
Name: close, dtype: float64
Didukung oleh

Ini menghasilkan Peringkat Volatilitas Bulanan Rata-Rata akhir kita. Secara numerik kita dapat melihat bahwa bulan ke-10 (Oktober) adalah yang tertinggi dan 12 (Desember) yang terendah.

#the final average results over 32 years
final
Didukung oleh

date
1     6.818182
2     6.666667
3     6.575758
4     7.303030
5     6.606061
6     6.030303
7     6.031250
8     5.875000
9     6.406250
10    7.531250
11    6.343750
12    5.218750
Name: close, dtype: float64
Didukung oleh

Memilih visualisasi yang tepat penting karena kita ingin menceritakan kisah data melalui bagan yang paling berdampak. Untuk analisis ini, kita ingin menampilkan dengan jelas volatilitas tertinggi dan terendah. Untuk itu, kita akan menggunakan diagram batang matplotlib dan memberi anotasi serta warna yang sesuai agar lebih berdampak.

 1. Dengan mengindeks b_plot (b_plot[9]) pada batang yang diinginkan, kita dapat mengatur warna untuk menyoroti nilai tertinggi dan terendah.

 2. Untuk menambahkan nilai AMVR ke batang masing-masing, kita enumerasi setiap nilai AMVR (variabel 'final') dan membulatkan hasil ke 2 angka desimal. 'i' dan 'v' mewakili 'index' dan 'value' untuk variabel 'final', jadi 'i' akan berkisar dari 1–12, dan 'value' memberi nilai AMVR dibulatkan 2 desimal. Kita gunakan variabel ini dalam plt.text(). Parameter pertama plt.text() adalah lokasi sumbu x untuk label, jadi kita teruskan 'i' untuk setiap batang dan offset 0,8 agar berada di tengah. Parameter kedua adalah versi string dari nilai 'final'.

 3. Untuk menampilkan nilai mean, kita gunakan axhline dan gunakan 'label' untuk menambahkan legenda di pojok kanan atas. Cukup teruskan nilai sumbu y, yaitu mean, lalu tambahkan gaya.

Ini memberi visual terbaik sejauh ini dari fenomena yang kita selidiki. Kita dapat melihat dengan jelas bahwa Oktober adalah bulan paling volatil dan Desember yang paling rendah. Penting, kita juga dapat melihat bahwa Desember adalah nilai yang lebih 'ekstrem' dalam istilah absolut. Ini akan penting untuk analisis kita di bagian berikutnya.

#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 in enumerate(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()
Didukung oleh

Itulah data kita, sekarang beralih ke pengujian Hipotesis…

Pengujian Hipotesis: Apa Pertanyaannya?

Pengujian hipotesis adalah salah satu teknik paling mendasar dalam ilmu data, namun juga salah satu yang paling menakutkan dan sering disalahpahami. Dasar dari ketakutan ini adalah cara pengajarannya di Stat 101, di mana kita diberi tahu untuk:

"lakukan uji t, apakah one-sided atau two-sided? Pilih statistik uji yang sesuai seperti uji t Welch, hitung derajat kebebasan, hitung skor t, lihat nilai kritis di tabel, bandingkan nilai kritis dengan statistik t …"

Wajar jika ini menimbulkan kebingungan tentang uji apa yang harus dilakukan dan bagaimana melakukannya. Namun, semua teknik statistik klasik untuk melakukan uji hipotesis dikembangkan pada saat daya komputasi sangat terbatas dan pada dasarnya merupakan solusi analitis bentuk tertutup untuk menghitung p-value, itu saja! Tetapi dengan komplikasi tambahan berupa keharusan memilih rumus yang tepat untuk situasi tertentu, karena asumsi mereka yang ketat dan terkadang buram.

Tapi bergembiralah!

Ada cara yang lebih baik. Simulasi.

Untuk memahami bagaimana simulasi dapat membantu, mari ingat kembali apa itu uji hipotesis:

Kita ingin menguji “apakah efek yang diamati pada data kita itu nyata atau apakah itu bisa terjadi hanya karena kebetulan” dan untuk melakukan uji ini kita melakukan hal berikut:

• Pilih ‘statistik uji’ yang sesuai: ini hanyalah angka yang mengukur efek yang diamati. Dalam kasus kita, kita akan memilih deviasi absolut AMVR dari mean.

• Menyusun Hipotesis Nol: ini hanyalah versi data di mana efek yang diamati tidak ada. Dalam kasus kita, kita akan mengacak label data berulang kali (permutasi). Justifikasinya dirinci di bawah.

• Hitung p-value: ini adalah probabilitas melihat efek yang diamati di antara data nol, dengan kata lain, karena kebetulan. Kita lakukan ini melalui simulasi berulang terhadap data nol. Dalam kasus kita, kita mengacak label ‘tanggal’ dari data berkali-kali dan cukup menghitung kejadian statistik uji saat muncul melalui banyak simulasi.

Itulah pengujian hipotesis dalam 3 langkah! Apa pun fenomena yang kita uji, pertanyaannya selalu sama: “apakah efek yang diamati itu nyata, atau karena kebetulan”

Hanya ada satu uji! Blog hebat oleh Allen Downey ini memuat lebih banyak detail tentang pengujian hipotesis

Kekuatan nyata simulasi adalah kita harus mengekspresikan secara eksplisit apa asumsi model kita melalui kode. Sedangkan teknik klasik bisa menjadi ‘kotak hitam’ dalam hal asumsi-asumsinya.

Contoh di bawah: Plot kiri menunjukkan data nyata dan efek yang diamati dengan probabilitas tertentu (hijau). Plot kanan adalah data nol simulasi kita dengan pencatatan kapan efek yang diamati terlihat karena kebetulan (merah). Inilah dasar pengujian hipotesis, berapa probabilitas melihat efek yang diamati pada data nol kita.

Bagian paling krusial dari pengujian hipotesis adalah memperjelas pertanyaan apa yang kita coba jawab. Dalam kasus kita kita bertanya:

“Mungkinkah nilai paling ekstrem terjadi karena kebetulan?”

Nilai paling ekstrem kita definisikan sebagai deviasi absolut AMVR terbesar dari mean. Pertanyaan ini membentuk hipotesis nol kita.

Dalam data kita, nilai paling ekstrem adalah nilai Desember (1,23) bukan nilai Oktober (1,08), karena kita melihat deviasi absolut paling signifikan dari mean, bukan sekadar volatilitas tertinggi.

 1. Untuk mendapatkan deviasi absolut, kita cukup mengurangkan setiap nilai dari mean dan menerapkan fungsi abs().

 2. Dengan menggunakan metode sort_values(), kita dapat mengurutkan hasil dari yang terkecil ke terbesar. Lalu pilih 2 nilai terbesar, Okt & Des.

#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)
Didukung oleh

date
9     0.044271
11    0.106771
3     0.125237
5     0.155540
2     0.216146
1     0.367661
7     0.419271
6     0.420218
8     0.575521
4     0.852509
10    1.080729
12    1.231771
Name: close, dtype: float64
Extreme Dec value: 1.231770833333333
Extreme Oct value: 1.080729166666667
Didukung oleh

Simulasi

Sekarang kita tahu pertanyaan apa yang diajukan, kita perlu menyusun ‘Model Nol’.

Ada beberapa opsi di sini:

• Model parametrik. Jika kita memiliki gambaran bagus tentang distribusi data atau sekadar membuat asumsi atasnya, kita bisa menggunakan teknik pengujian hipotesis ‘klasik’, t-test, X², one-way ANOVA, dll. Model-model ini bisa restriktif dan agak seperti kotak hitam jika peneliti tidak sepenuhnya memahami asumsinya.

• Simulasi langsung. Kita dapat membuat asumsi tentang proses pembangkitan data dan mensimulasikannya langsung. Misalnya, kita bisa menetapkan model deret waktu ARMA untuk data keuangan yang kita modelkan dan dengan sengaja merancangnya agar tidak memiliki musiman. Ini bisa menjadi pilihan yang masuk akal untuk masalah kita. Namun, jika kita tahu proses pembangkitan data untuk S&P500, kita pasti sudah kaya!

• Simulasi melalui Resampling. Ini pendekatan yang akan kita ambil. Dengan berulang kali mengambil sampel acak dari dataset yang ada dan mengacak label, kita dapat membuat efek yang diamati sama kemungkinannya di antara semua label dalam data kita (dalam kasus kita label adalah tanggal), sehingga memberikan dataset nol yang diinginkan.

Sampling adalah topik besar, tetapi kita akan fokus pada satu teknik khusus, permutasi atau pengacakan.

Untuk mendapatkan model nol yang diinginkan, kita perlu menyusun dataset yang tanpa musiman. Jika hipotesis nol benar, bahwa tidak ada musiman dalam data dan efek yang diamati terjadi hanya karena kebetulan, maka label untuk setiap bulan (Jan, Feb, dll.) tidak bermakna dan karena itu kita bisa mengacak data berulang kali untuk membangun apa yang oleh statistik klasik disebut ‘sebaran sampling dari statistik uji di bawah hipotesis nol’. Ini memiliki efek yang diinginkan yaitu membuat fenomena yang diamati (nilai Desember ekstrem) sama kemungkinannya untuk semua bulan, yang persis seperti yang dibutuhkan model nol kita.

Untuk membuktikan betapa kuatnya teknik simulasi dengan daya komputasi modern, kode dalam contoh ini sebenarnya akan melakukan permutasi pada data harga harian, yang membutuhkan daya proses jauh lebih banyak, namun tetap selesai dalam hitungan detik pada CPU modern.

Catatan: mengacak label harian atau bulanan sama-sama akan memberi dataset nol yang diinginkan dalam kasus kita.

                          Shuffle 'date' label to create null dataset
Didukung oleh

Sumber belajar hebat tentang sampling oleh Julian Simon.

Catatan: Cara uji kita dibangun setara dengan uji dua sisi menggunakan metode ‘klasik’, seperti uji t Welch atau ANOVA, dll., karena kita tertarik pada nilai paling ekstrem, baik di atas maupun di bawah mean.

Keputusan-keputusan ini adalah pilihan desain, dan kita memiliki kebebasan ini karena Model Nol hanyalah itu, sebuah Model! Artinya kita bisa menentukan parameternya sesuka kita, kuncinya adalah benar-benar jelas pertanyaan apa yang ingin kita jawab.

Jadi kita ingin menggunakan Python untuk mensimulasikan banyak data guna membuat dataset nol. Untuk mencapainya, kita akan mensimulasikan 1000 set berisi 12 AMVR, mengacak label 'tanggal' setiap kali untuk membangun sebaran sampling. Keluaran dari kode ini disertakan di bawah dalam bagian p-hacking.

 1. Pertama kita mendefinisikan beberapa penampung untuk menyimpan hasil: kita akan menggunakan pd.Dataframe() dan array sederhana [ ].

 2. Selanjutnya, kita mendefinisikan pencacah dan setel ke nol, lalu memulai for-loop 1000 iterasi untuk membuat data simulasi.

 3. Baris pertama di dalam for-loop mengambil imbal hasil harian asli (dari awal tutorial) dan menggunakan metode pandas sample(). Ini mengambil sampel acak dari data imbal hasil harian sejumlah kali tertentu. Dalam kasus kita, ada 8191 titik data (252 hari perdagangan selama 32 tahun). Kita juga membuang indeks menggunakan reset_index(), sehingga kita dapat menambahkan indeks baru. Kita perlu melakukan ini karena indeks tanggal asli ikut teracak bersama data, sedangkan kita ingin 'mengacak' hanya datanya.

 4. Baris berikutnya menambahkan indeks tanggal baru dengan menetapkan ulang indeks daily_ret_shuffle dengan pd.bdate_range yang panjangnya sama dengan data asli. Kita gunakan bdate_range alih-alih date_range, karena kita ingin hari kerja (5, Sen - Jum) bukan 7 hari per minggu.

 5. Sekarang setelah kita menghasilkan data 'teracak', kita melakukan pengolahan data yang sama seperti di awal untuk menyusun nilai AMVR. Itulah yang dilakukan 3 baris kode berikutnya.

 6. Setelah memiliki nilai AMVR simulasi, kita perlu menambahkannya ke dataframe yang akan menyimpan semua 1000 kali simulasi. pd.concat mengambil dataframe saat ini dan menambahkan setiap dataframe berikutnya di bagian akhir. Kita pilih axis 1, sehingga kita mendapatkan kolom data.

 7. Maximonth adalah variabel yang menyimpan hanya nilai tertinggi dari setiap simulasi (akan kita gunakan nanti dalam analisis).

 8. Karena analisis kita memerlukan nilai AMVR absolut, 3 baris kode berikutnya mengambil seluruh 1000 kali AMVR dan meratakannya menjadi satu array besar. Lalu kita ambil mean dan kurangkan nilai individual dari mean ini dan akhirnya ambil nilai absolut.

 9. Kita juga melakukan hal yang sama untuk data 'tertinggi saja'. Catatan: di sini kita menggunakan list alih-alih dataframe yang berarti kita bisa menggunakan list comprehension untuk menghitung abs(AMVR) untuk masing-masing nilai tertinggi kita.

Sekarang kita memiliki semua data yang dibutuhkan untuk menyelesaikan analisis, baik observasi asli maupun dataset simulasi baru kita.

#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=1000
for i in range(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]
Didukung oleh

p-Hacking

Bagian menariknya di sini. Kita telah menyusun hipotesis untuk diuji, kita telah menghasilkan data simulasi dengan mengacak label ‘tanggal’ pada data, sekarang kita perlu melakukan uji hipotesis untuk mencari probabilitas mengamati hasil yang setinggi hasil Desember dengan asumsi hipotesis nol (tanpa musiman) benar.

Sebelum melakukan uji, mari tetapkan ekspektasi kita.

Berapa probabilitas melihat setidaknya satu hasil signifikan dengan tingkat signifikansi 5%?

= 1-p(tidak signifikan)

= 1-(1–0.05)¹²

= 0.46

jadi ada peluang 46% untuk melihat setidaknya satu bulan dengan hasil signifikan, dengan asumsi hipotesis nol kita benar.

Sekarang mari bertanya, untuk setiap uji individual (membandingkan AMVR absolut 12 bulan dengan mean) berapa banyak nilai signifikan yang seharusnya kita harapkan muncul di antara data acak tanpa musiman?

12 x 0,05 = 0,6

Jadi dengan tingkat signifikansi 0,05, kita seharusnya mengharapkan tingkat positif palsu sebesar 0,6. Dengan kata lain, untuk setiap uji (dengan data nol) yang membandingkan semua 12 AMVR bulan terhadap mean, 0,6 bulan akan menunjukkan hasil signifikan. (jelas, kita tidak bisa memiliki kurang dari 1 bulan yang menunjukkan hasil, tetapi dalam pengujian berulang matematika akan menuju angka ini).

Sepanjang pekerjaan ini, kita menekankan pentingnya benar-benar jelas dengan pertanyaan yang coba kita jawab. Masalah dengan ekspektasi yang baru saja kita hitung adalah kita berasumsi kita menguji hasil signifikan terhadap semua 12 bulan! Itulah mengapa probabilitas melihat setidaknya satu positif palsu begitu tinggi di 46%.

Ini adalah contoh bias perbandingan ganda di mana kita memperluas ruang pencarian dan meningkatkan kemungkinan menemukan hasil signifikan. Ini menjadi masalah karena kita bisa menyalahgunakan fenomena ini untuk memilih parameter model yang memberi p-value yang ‘diinginkan’.

Inilah esensi p-Hacking

Untuk menggambarkan dampak p-hacking dan cara mengurangi multiplicity, kita perlu memahami perbedaan halus tetapi signifikan antara 2 pertanyaan berikut:

• “Berapa probabilitas Desember akan tampak se-ekstrem ini karena kebetulan?”
• “Berapa probabilitas bulan apa pun akan tampak se-ekstrem ini karena kebetulan?”

Keindahan simulasi terletak pada kesederhanaannya. Kode berikut adalah semua yang kita butuhkan untuk menghitung p-value guna menjawab pertanyaan pertama. Kita cukup menghitung berapa banyak nilai dalam dataset kita menggunakan semua 12000 deviasi AMVR (12 bulan x 1000 percobaan) yang lebih besar daripada nilai Desember yang diamati. Kita mendapatkan p-value 4,4%, mendekati batas 5% yang agak sewenang-wenang, tetapi tetap signifikan.

#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=0
for i in abs_all_months:
    if i> Dec_value:
        count+=1
ans = count/len(abs_all_months)        
print('p-value:', ans )
Didukung oleh

p-value: 0.04425
Didukung oleh

Untuk menjawab pertanyaan kedua dan menghindari multiplicity, alih-alih membandingkan hasil kita dengan distribusi yang dibuat dari semua 12000 deviasi AMVR, kita hanya mempertimbangkan nilai tertinggi dari masing-masing 1000 percobaan AMVR absolut. Ini menghasilkan p-value 23%, jelas tidak signifikan!

#same again but just considering highest AMVR for each of 100 trials
count=0
for i in abs_highest:
    if i> Dec_value:
        count+=1
ans = count/len(abs_highest)        
print('p-value:', ans )
Didukung oleh

p-value: 0.236
Didukung oleh

Sekarang setelah kita memiliki hasil akhir, mari plot distribusi ini untuk menunjukkan efek p-Hacking dan hasil analisis kita secara visual:

 1. Pertama kita gunakan np.quantile() untuk mencari tingkat signifikansi 5%, kita plot ini menggunakan axvline pada plot bawah.

 2. Selanjutnya kita definisikan 4 subplot kotak. Di sini plt.subplots mengembalikan 'fig' = figur kita dan 'ax1,ax2,ax3,ax4' = variabel yang merepresentasikan masing-masing dari 4 subplot.

 3. #plot 1 menampilkan kolom pertama. Di sini kita cukup mendefinisikan histogram data 'abs_all_months' dan pilih type='bar'. Untuk plot bawah kita set cumulative='True', ini memberi kita cdf alih-alih pdf. Bins menentukan berapa banyak batang yang kita miliki untuk setiap nilai yang diisikan, 30 adalah angka yang masuk akal untuk data kita. Lalu kita format plot menggunakan teknik yang sudah kita pelajari, seperti axvline untuk memplot nilai signifikansi dan observasi.

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)
Didukung oleh

<matplotlib.legend.Legend at 0x280b4eb0b00>
Didukung oleh

Kolom kiri adalah data yang menjawab pertanyaan 1 dan kolom kanan, pertanyaan 2. Baris atas adalah distribusi probabilitas, dan baris bawah adalah CDF. Garis putus-putus merah adalah tingkat signifikansi 5% yang kita tetapkan secara sewenang-wenang. Garis biru adalah nilai AMVR Desember ekstrem asli sebesar 1,23.

Plot sisi kiri menunjukkan bahwa nilai Desember asli signifikan pada level 5%, tetapi nyaris! Namun, saat kita memperhitungkan bias perbandingan ganda, pada plot kanan ambang batas signifikansi naik dari sekitar 1,2 (abs AMVR) ke sekitar 1,6 (lihat garis merah).

Dengan memperhitungkan bias perbandingan ganda, nilai Desember 1,23 kita tidak lagi signifikan!

Dengan mempertimbangkan pertanyaan spesifik yang kita coba jawab dan menghindari bias perbandingan ganda, kita menghindari p-hacking pada model dan menghindari menampilkan hasil signifikan saat sebenarnya tidak ada.

Untuk mengeksplorasi lebih jauh p-hacking dan bagaimana ia dapat disalahgunakan untuk menceritakan kisah tertentu tentang data kita, lihat aplikasi interaktif hebat dari FiveThirtyEight

Kesimpulan

• Kita telah mempelajari bahwa pengujian hipotesis bukanlah momok besar yang kita kira. Cukup ikuti 3 langkah di atas untuk menyusun model Anda untuk jenis data atau statistik uji apa pun.

• Kita telah menunjukkan bahwa mengajukan pertanyaan yang tepat sangat penting untuk analisis ilmiah. Perubahan kecil dalam perumusan dapat menghasilkan model yang sangat berbeda dengan hasil yang sangat berbeda.

• Mudah-mudahan, kekuatan fungsionalitas Python tingkat menengah dan kemampuannya memberdayakan Anda melakukan pengujian statistik kini terlihat jelas! Hanya dengan beberapa baris kode, kita telah menyusun dan menguji fenomena dunia nyata dan mampu menarik kesimpulan yang dapat ditindaklanjuti.

• Kita membahas pentingnya mengenali dan mengoreksi bias perbandingan ganda serta menghindari jebakan p-hacking dan menunjukkan bagaimana hasil yang tampaknya signifikan dapat menjadi tidak signifikan.

• Dengan semakin banyaknya ‘big data’ bersama tekanan akademis untuk menghasilkan makalah riset dengan temuan ‘baru’ atau tekanan politis untuk menunjukkan hasil sebagai ‘signifikan’, godaan untuk melakukan p-hacking semakin besar. Dengan belajar mengenali kapan kita melakukannya dan mengoreksinya, kita dapat menjadi peneliti yang lebih baik dan pada akhirnya menghasilkan hasil ilmiah yang lebih akurat dan karenanya dapat ditindaklanjuti!

Catatan Penulis: Hasil kami sedikit berbeda dari riset Winton asli, sebagian karena memiliki dataset yang agak berbeda (32 tahun vs. 87 tahun) dan mereka menjadikan Oktober sebagai bulan yang menarik sedangkan kami Desember. Selain itu, mereka menggunakan metode yang tidak diungkapkan untuk ‘data simulasi’ mereka sementara kami telah menjelaskan secara eksplisit, melalui kode, metodologi untuk membuat data tersebut. Kami telah membuat sejumlah asumsi pemodelan sepanjang pekerjaan ini, lagi-lagi, ini dibuat eksplisit dan dapat dilihat dalam kode. Pilihan desain dan pemodelan ini merupakan bagian dari proses ilmiah, selama dibuat eksplisit, analisis tetap bernilai.

Jika Anda ingin belajar lebih banyak tentang keuangan di Python, ikuti kursus DataCamp Intro to Python for Finance dan lihat Python For Finance Tutorial: Algorithmic Trading.

Ikuti saya di twitter.com/pdquant untuk info terbaru!