Tìm hiểu cách xử lý dữ liệu chuỗi thời gian với pandas và thực hiện kiểm định ý nghĩa thống kê bằng mô phỏng trong Python, để phân tích biến động thị trường chứng khoán.
Tháng 10 trong lịch sử là tháng biến động mạnh nhất đối với cổ phiếu, nhưng đây là tín hiệu bền vững hay chỉ là nhiễu trong dữ liệu?
"Trong 32 năm qua, tháng 10 là tháng biến động nhất tính trung bình đối với S&P500 và tháng 12 là tháng ít biến động nhất".
Trong hướng dẫn này, chúng ta sẽ dùng Python để đi qua toàn bộ quá trình phân tích và kiểm định hiện tượng này nhằm xác định liệu nó có ý nghĩa thống kê hay không.
Chúng ta sẽ sử dụng các chức năng chuỗi thời gian của pandas, như Resample, để biến đổi dữ liệu giá cổ phiếu thô thành định dạng làm dữ liệu trở nên trực quan.
Sau đó, chúng ta sẽ tìm hiểu cách thực hiện kiểm định giả thuyết bằng mô phỏng trong Python thay vì dùng các công thức rườm rà.
Cuối cùng, chúng ta sẽ bàn về một vấn đề lớn trong phân tích thống kê, hiện tượng sai lệch so sánh nhiều lần (Multiple Comparison Bias), cách xử lý và minh họa trực quan qua matplotlib tác động của 'p-Hacking'.
Mục tiêu của chúng ta:
Minh họa cách dùng Pandas để phân tích chuỗi thời gian
Hiểu cách xây dựng một kiểm định giả thuyết
Dùng mô phỏng với Python để thực hiện kiểm định giả thuyết
Cho thấy tầm quan trọng của việc tính đến sai lệch so sánh nhiều lần
Dữ liệu của chúng ta:
Chúng ta sẽ sử dụng dữ liệu S&P500 theo ngày cho phân tích này, cụ thể là giá đóng cửa thô hằng ngày từ 1986 đến 2018 (khá khó tìm, nên tôi đã công khai).
Cảm hứng cho bài viết đến từ Winton, chúng ta sẽ tái hiện lại ở đây, dù chỉ với 32 năm dữ liệu so với 87 năm của họ.
Xử lý dữ liệu với Pandas
Để trả lời câu hỏi liệu biến động cực đoan ở một số tháng có thực sự đáng kể và do đó có khả năng tiếp diễn hay không, chúng ta cần biến đổi 32 năm dữ liệu giá thành một định dạng thể hiện rõ hiện tượng đang khảo sát.
Định dạng lựa chọn là bảng xếp hạng biến động trung bình theo tháng (AMVR).
Đoạn mã sau cho thấy cách đưa dữ liệu giá thô về định dạng này. Bắt đầu thôi!
Đầu tiên là các import tiêu chuẩn. (matplotlib.patches cho phép chúng ta kiểm soát việc tạo kiểu cho từng cột trong biểu đồ cột)
#standard importsimport numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
%matplotlib inline
Được hỗ trợ bởi
Một mẹo hay để biểu đồ hiển thị toàn chiều rộng trong Jupyter notebooks:
#resize charts to fit screen if using Jupyter Notebook
plt.rcParams['figure.figsize']=[15,5]Được hỗ trợ bởi
Import dữ liệu bằng phương thức read_csv nhận một đường dẫn, ở đây là url đến dữ liệu. Ta cũng chỉ định dùng 'date' làm chỉ mục và tự động phân tích cú pháp ngày tháng từ văn bản cho phù hợp nhất có thể.
#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()Được hỗ trợ bởi
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
Điều này cho ta giá đóng cửa thô chưa điều chỉnh của S&P500 (SPX). Giờ ta cần chuyển giá thô sang lợi nhuận % theo ngày. Có hai lựa chọn:
Lấy log tự nhiên của giá. Điều này sẽ cho xấp xỉ lợi nhuận ngày thực.
Dùng phương thức 'pct_change()' của pandas để tính trực tiếp phần trăm thay đổi hằng ngày.
Với mục đích của chúng ta, sẽ dùng cách 2, vì pandas có thể xử lý tức thì tập dữ liệu cỡ này (hơn 8000 giá trị). Ta cũng sẽ làm sạch dữ liệu bằng cách bỏ giá trị đầu tiên trở thành 'NaN' do không có thay đổi giá so với ngày trước đó. pct_change() nhận tham số tùy chọn 'periods' để đổi khoảng dịch; ta giữ mặc định là 1.
#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()Được hỗ trợ bởi
date
1986-01-030.0061551986-01-06-0.0010911986-01-070.0149541986-01-08-0.0272681986-01-09-0.008944Name: close,dtype: float64
Được hỗ trợ bởi
Bước tiếp theo là lấy các giá trị % thay đổi theo ngày và chuyển chúng thành 'biến động hàng tháng đã quy đổi theo năm'. Dòng đầu tiên dưới đây cho thấy ta có thể thực hiện chuyển đổi này chỉ bằng một dòng mã. Đó là sức mạnh của Pandas! Nhưng hãy bóc tách từng bước.
1. Để có mnthly_annu, đầu tiên ta dùng phương thức 'resample' trên lợi nhuận ngày. Resample cho phép thay đổi tần suất kỳ dữ liệu. Nó nhận một 'chuỗi ký hiệu tần suất', nói đơn giản là một chữ cái tương ứng tần suất mong muốn. Bao gồm:
B tần suất ngày làm việc
C tần suất ngày làm việc tùy chỉnh
D tần suất ngày lịch
W tần suất tuần
M tần suất cuối tháng
Vì muốn lấy mẫu lại từ lợi nhuận ngày sang tháng, ta truyền 'M' làm tham số.
2. Tiếp theo cần quyết định cách tính giá trị hàng tháng mới, ví dụ có thể cộng dồn, nhân, v.v. Với phân tích của ta, cần thước đo biến động và độ lệch chuẩn phù hợp, nên ta nối thêm std() vào dữ liệu đã resample để lấy biến động tháng.
3. Bước cuối là quy đổi theo năm con số này. Ta làm đơn giản bằng cách nhân với căn bậc hai của 12, với 12 là số kỳ (tháng) trong năm. Điều này cho ta các giá trị biến động tháng quy đổi theo năm cần thiết.
Thực hành tốt trong phân tích dữ liệu là trực quan hóa dữ liệu khi tiến hành phân tích. Hãy xem dữ liệu biến động đã quy đổi theo năm của chúng ta.
Biểu đồ dưới cho thấy rõ các sự kiện lớn trên thị trường, như Black Monday và Khủng hoảng Tài chính 2008. Phương thức 'axvspan' của matplotlib cho phép thêm các vạch dọc (axhspan là phương thức tương ứng cho phương ngang), nhận 'xmin' và 'xmax' làm tham số, xác định bề rộng của hình chữ nhật tạo ra; vì trục x là chỉ mục datetime, ta có thể truyền trực tiếp các năm muốn làm nổi bật. Tham số alpha điều chỉnh độ trong suốt của khối màu để vẫn thấy biểu đồ bên dưới.
Thuộc tính mpatches cho phép tạo một hộp chú giải tùy chỉnh. Ta làm bằng cách định nghĩa biến 'labs' với màu, alpha và văn bản mong muốn, sau đó truyền biến này vào tham số 'handles' trong plt.legend để hiển thị chú giải.
#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])Được hỗ trợ bởi
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>Được hỗ trợ bởi
Vậy là chúng ta đã thấy một trong những phương thức mạnh mẽ của pandas, resample. Giờ hãy dùng một phương thức khác, 'groupby'. Chúng ta cần đi từ các giá trị biến động tháng đã quy đổi theo năm sang thước đo AMVR mong muốn. Lại một lần nữa, ta có thể đạt được chỉ với vài dòng mã.
1. Đầu tiên, áp dụng 'groupby' lên mnthly_annu. Groupby cần một tham số xác định cách nhóm dữ liệu. Có thể là một hàm hoặc như ở đây, một Series. Ta truyền 'mnthly_annu.index.year', đơn giản là thuộc tính năm trên chỉ mục datetime của mnthly_annu. Điều này sẽ nhóm các giá trị biến động tháng cho từng năm trong bộ dữ liệu 32 năm.
2. Tiếp theo, áp dụng phương thức 'rank', sắp xếp dữ liệu theo thứ tự tăng dần.
3. Cuối cùng, lặp lại và lấy trung bình theo tất cả các năm cho từng tháng. Điều này cho ra các giá trị AMVR cuối cùng!
#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()Được hỗ trợ bởi
count 12.000000
mean 6.450521
std 0.627458
min 5.21875025%6.03101350%6.49100475%6.704545
max 7.531250Name: close,dtype: float64
Được hỗ trợ bởi
Điều này cho ta Bảng xếp hạng Biến động Trung bình theo Tháng cuối cùng. Về số liệu, ta thấy tháng 10 là cao nhất và tháng 12 là thấp nhất.
#the final average results over 32 years
final
Được hỗ trợ bởi
date
16.81818226.66666736.57575847.30303056.60606166.03030376.03125085.87500096.406250107.531250116.343750125.218750Name: close,dtype: float64
Được hỗ trợ bởi
Chọn trực quan hóa phù hợp cho dữ liệu là quan trọng vì chúng ta muốn kể câu chuyện dữ liệu qua những biểu đồ tác động nhất. Với phân tích này, ta muốn thể hiện rõ mức biến động cao nhất và thấp nhất. Ta sẽ dùng biểu đồ cột của matplotlib và chú thích, tô màu phù hợp để tăng tác động.
1. Bằng cách lập chỉ mục b_plot (b_plot[9]) với cột mong muốn, ta có thể đặt màu để làm nổi bật giá trị cao nhất và thấp nhất.
2. Để thêm giá trị AMVR lên cột tương ứng, ta duyệt liệt kê qua từng giá trị AMVR (biến 'final') và làm tròn đến 2 chữ số thập phân. 'i' và 'v' lần lượt là 'index' và 'value' của biến 'final', nên 'i' sẽ từ 1-12, và 'value' là giá trị AMVR làm tròn 2 chữ số. Ta dùng các biến này trong plt.text(). Tham số thứ nhất của plt.text() là vị trí trục x cho nhãn, nên ta truyền 'i' cho mỗi cột và dịch 0.8 để căn giữa. Tham số thứ hai là chuỗi biểu diễn giá trị 'final'.
3. Để hiển thị giá trị trung bình, ta dùng axhline và dùng 'label' để thêm chú giải ở góc trên bên phải. Chỉ cần truyền giá trị trục y, tức giá trị trung bình, rồi thêm định dạng.
Biểu đồ này là trực quan tốt nhất cho hiện tượng đang khảo sát. Ta thấy rõ tháng 10 biến động mạnh nhất và tháng 12 thấp nhất. Quan trọng là ta cũng thấy tháng 12 là giá trị 'cực trị' hơn, theo giá trị tuyệt đối. Điều này sẽ quan trọng cho phân tích ở phần tiếp theo.
#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()Được hỗ trợ bởi
Vậy là xong phần dữ liệu, giờ đến Kiểm định giả thuyết…
Kiểm định giả thuyết: Câu hỏi là gì?
Kiểm định giả thuyết là một trong những kỹ thuật nền tảng của khoa học dữ liệu, nhưng cũng là một trong những thứ gây e ngại và dễ bị hiểu sai. Căn nguyên nằm ở cách dạy trong môn Thống kê cơ bản, nơi ta được bảo:
"thực hiện t-test, đó là kiểm định một phía hay hai phía? Chọn thống kê kiểm định phù hợp như t-test của Welch, tính bậc tự do, tính điểm t, tra giá trị tới hạn trong bảng, so sánh giá trị tới hạn với thống kê t ……"
Dễ hiểu khi điều này dẫn đến bối rối về việc nên dùng kiểm định nào và tiến hành ra sao. Tuy nhiên, tất cả các kỹ thuật thống kê cổ điển để thực hiện kiểm định giả thuyết được phát triển vào thời điểm khả năng tính toán rất hạn chế và chỉ là các nghiệm đóng để tính p-value, hết! Nhưng kèm thêm sự phức tạp là phải chọn đúng công thức cho tình huống cụ thể, do các giả định hạn chế và đôi khi mù mờ.
Nhưng hãy mừng lên!
Có cách tốt hơn. Mô phỏng.
Để hiểu mô phỏng giúp gì, hãy nhắc lại kiểm định giả thuyết là gì:
Chúng ta muốn kiểm tra “liệu hiệu ứng quan sát được trong dữ liệu là thật hay có thể xảy ra đơn thuần do ngẫu nhiên” và để thực hiện, ta làm như sau:
Chọn một ‘thống kê kiểm định’ phù hợp: đơn giản là một con số đo lường hiệu ứng quan sát được. Trong trường hợp này, ta chọn độ lệch tuyệt đối của AMVR so với trung bình.
Xây dựng Giả thuyết Không: đây chỉ là một phiên bản dữ liệu mà ở đó hiệu ứng quan sát không hiện diện. Ở đây, ta sẽ xáo trộn nhãn của dữ liệu lặp đi lặp lại (permutation). Lý do được trình bày bên dưới.
Tính p-value: đây là xác suất thấy hiệu ứng quan sát được trong dữ liệu null, nói cách khác là do ngẫu nhiên. Ta làm điều này thông qua mô phỏng lặp lại dữ liệu null. Ở đây, ta xáo trộn nhãn ‘date’ của dữ liệu nhiều lần và đơn giản đếm số lần thống kê kiểm định xuất hiện qua nhiều lần mô phỏng.
Đó là kiểm định giả thuyết trong 3 bước! Bất kể ta kiểm tra hiện tượng nào, câu hỏi luôn như nhau: “hiệu ứng quan sát được là thật, hay do ngẫu nhiên”
Sức mạnh thực sự của mô phỏng là buộc chúng ta phải diễn đạt rõ ràng các giả định mô hình qua mã. Trong khi kỹ thuật cổ điển có thể là một ‘hộp đen’ về các giả định của chúng.
Ví dụ dưới đây: Biểu đồ trái là dữ liệu thật và hiệu ứng quan sát với một xác suất nhất định (màu xanh lá). Biểu đồ phải là dữ liệu null mô phỏng với ghi nhận khi nào hiệu ứng quan sát xuất hiện do ngẫu nhiên (màu đỏ). Đây là nền tảng của kiểm định giả thuyết: xác suất thấy hiệu ứng quan sát trên dữ liệu null là bao nhiêu.
Phần quan trọng nhất của kiểm định giả thuyết là rõ ràng về câu hỏi cần trả lời. Ở đây chúng ta hỏi:
“Giá trị cực trị nhất có thể xảy ra do ngẫu nhiên không?”
Giá trị cực trị nhất được định nghĩa là độ lệch tuyệt đối của AMVR so với trung bình lớn nhất. Câu hỏi này tạo nên giả thuyết không của chúng ta.
Trong dữ liệu, giá trị cực trị nhất là giá trị tháng 12 (1.23) chứ không phải tháng 10 (1.08), vì ta xét độ lệch tuyệt đối lớn nhất so với trung bình, không chỉ là biến động cao nhất.
1. Để có độ lệch tuyệt đối, chỉ cần lấy mỗi giá trị trừ đi trung bình và áp dụng hàm abs().
2. Dùng sort_values() để sắp xếp kết quả từ nhỏ đến lớn. Sau đó chọn ra 2 giá trị lớn nhất, Tháng 10 & Tháng 12.
#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)Được hỗ trợ bởi
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.080729166666667Được hỗ trợ bởi
Mô phỏng
Giờ ta đã biết câu hỏi cần hỏi, cần xây dựng ‘Mô hình Null’.
Có một số lựa chọn:
Mô hình tham số. Nếu ta có hiểu biết tốt về phân phối dữ liệu hoặc đơn giản đặt giả định về nó, có thể dùng các kỹ thuật kiểm định giả thuyết ‘cổ điển’, t-test, X², ANOVA một nhân tố, v.v. Các mô hình này có thể hạn chế và phần nào là hộp đen nếu nhà nghiên cứu không hiểu rõ các giả định.
Mô phỏng trực tiếp. Ta có thể đặt giả định về quá trình sinh dữ liệu và mô phỏng trực tiếp. Ví dụ, chỉ định một mô hình chuỗi thời gian ARMA cho dữ liệu tài chính đang mô hình hóa và cố ý thiết kế để không có tính mùa vụ. Đây có thể là lựa chọn hợp lý cho vấn đề của ta. Tuy nhiên, nếu biết quá trình sinh dữ liệu của S&P500, có lẽ ta đã giàu rồi!
Mô phỏng qua lấy mẫu lại. Đây là cách ta sẽ dùng. Bằng cách lặp lại việc lấy mẫu ngẫu nhiên từ tập dữ liệu hiện có và xáo trộn nhãn, ta có thể làm cho hiệu ứng quan sát được có khả năng như nhau giữa mọi nhãn trong dữ liệu (ở đây nhãn là ngày), từ đó tạo ra tập dữ liệu null mong muốn.
Lấy mẫu là một chủ đề lớn, nhưng ta sẽ tập trung vào một kỹ thuật, hoán vị hay xáo trộn.
Để có mô hình null mong muốn, ta cần dựng một tập dữ liệu không có tính mùa vụ. Nếu giả thuyết không đúng, rằng không có mùa vụ trong dữ liệu và hiệu ứng quan sát chỉ do ngẫu nhiên, thì nhãn cho từng tháng (Jan, Feb, v.v.) là vô nghĩa và do đó ta có thể xáo trộn dữ liệu lặp lại để xây dựng thứ mà thống kê cổ điển gọi là ‘phân phối lấy mẫu của thống kê kiểm định dưới giả thuyết không’. Điều này khiến hiện tượng quan sát (giá trị tháng 12 cực trị) có khả năng như nhau cho mọi tháng, đúng như yêu cầu của mô hình null.
Để chứng minh sức mạnh của kỹ thuật mô phỏng với năng lực tính toán hiện đại, mã trong ví dụ này thực sự hoán vị dữ liệu giá theo ngày, đòi hỏi nhiều sức mạnh xử lý hơn, nhưng vẫn hoàn tất trong vài giây trên CPU hiện đại.
Lưu ý: xáo trộn nhãn theo ngày hoặc theo tháng đều cho tập dữ liệu null mong muốn trong trường hợp này.
Một tài nguyên tuyệt vời để học về lấy mẫu đến từ Julian Simon.
Lưu ý: Cách ta xây dựng kiểm định tương đương kiểm định hai phía theo phương pháp ‘cổ điển’, như t-test của Welch hay ANOVA, v.v., vì ta quan tâm đến giá trị cực trị nhất, ở trên hay dưới trung bình.
Những quyết định này là lựa chọn thiết kế, và ta có quyền tự do đó vì Mô hình Null chỉ là mô hình! Điều này nghĩa là ta có thể chỉ định tham số theo ý, điều mấu chốt là phải thật rõ ràng câu hỏi cần trả lời.
Vậy ta muốn dùng Python để mô phỏng nhiều dữ liệu nhằm tạo tập dữ liệu null. Để làm điều này, ta sẽ mô phỏng 1000 bộ gồm 12 AMVR, hoán vị nhãn 'date' mỗi lần để xây dựng phân phối lấy mẫu. Kết quả từ mã này được đưa ở phần p-hacking bên dưới.
1. Đầu tiên ta định nghĩa vài cấu trúc chứa kết quả: sẽ dùng một pd.Dataframe() và một mảng đơn [ ].
2. Tiếp theo, định nghĩa một biến đếm và đặt 0, rồi bắt đầu vòng lặp for 1000 lần để tạo dữ liệu mô phỏng.
3. Dòng đầu trong vòng lặp lấy lợi nhuận ngày gốc (từ đầu hướng dẫn) và dùng phương thức sample() của pandas. Phương thức này lấy mẫu ngẫu nhiên từ dữ liệu lợi nhuận ngày một số lần nhất định. Ở đây, ta có 8191 điểm dữ liệu (252 ngày giao dịch trong 32 năm). Ta cũng bỏ chỉ mục bằng reset_index(), để có thể thêm chỉ mục mới. Cần làm vậy vì chỉ mục ngày gốc sẽ bị xáo trộn cùng dữ liệu, trong khi ta muốn 'xáo trộn' chỉ dữ liệu.
4. Dòng tiếp theo thêm chỉ mục ngày mới bằng cách gán lại chỉ mục của daily_ret_shuffle với một pd.bdate_range có độ dài bằng dữ liệu gốc. Ta dùng bdate_range thay vì date_range vì cần ngày làm việc (5, Thứ Hai - Thứ Sáu) thay vì 7 ngày một tuần.
5. Giờ đã tạo dữ liệu 'xáo trộn', ta thực hiện lại quá trình xử lý như lúc đầu để tạo giá trị AMVR. Ba dòng mã tiếp theo làm việc này.
6. Khi đã có giá trị AMVR mô phỏng, ta cần nối chúng vào dataframe lưu trữ toàn bộ 1000 lượt mô phỏng. pd.concat nhận dataframe hiện tại và nối từng dataframe tiếp theo vào cuối. Ta chọn axis 1, để có các cột dữ liệu.
7. Maximonth là biến chỉ lưu giá trị cao nhất từ mỗi lần mô phỏng (sẽ dùng sau trong phân tích).
8. Vì phân tích cần giá trị AMVR tuyệt đối, 3 dòng mã tiếp theo lấy toàn bộ 1000 lượt AMVR và làm phẳng thành một mảng lớn. Sau đó lấy trung bình và trừ từng giá trị cho trung bình này rồi lấy giá trị tuyệt đối.
9. Ta cũng làm tương tự cho dữ liệu 'chỉ cao nhất'. Lưu ý: ở đây dùng list thay vì dataframe nên có thể dùng list comprehension để tính abs(AMVR) cho từng giá trị cao nhất.
Giờ ta đã có mọi dữ liệu cần để hoàn tất phân tích, cả quan sát gốc lẫn tập dữ liệu mô phỏng mới.
#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]Được hỗ trợ bởi
p-Hacking
Đây là phần thú vị. Ta đã xây dựng giả thuyết để kiểm định, đã tạo dữ liệu mô phỏng bằng cách xáo trộn nhãn ‘date’, giờ cần thực hiện kiểm định để tìm xác suất quan sát được kết quả đáng kể như kết quả tháng 12 với giả thuyết không (không có mùa vụ) là đúng.
Trước khi kiểm định, hãy đặt kỳ vọng.
Xác suất thấy ít nhất một kết quả có ý nghĩa với mức ý nghĩa 5% là bao nhiêu?
= 1-p(không có ý nghĩa)
= 1-(1–0.05)¹²
= 0.46
vậy có 46% khả năng thấy ít nhất một tháng có kết quả có ý nghĩa, giả thuyết không là đúng.
Giờ hãy hỏi, với từng kiểm định riêng lẻ (so sánh độ lệch tuyệt đối AMVR của 12 tháng so với trung bình) chúng ta kỳ vọng thấy bao nhiêu giá trị có ý nghĩa trong dữ liệu ngẫu nhiên, không mùa vụ?
12 x 0.05 = 0.6
Vậy với mức ý nghĩa 0.05, ta kỳ vọng tỷ lệ dương tính giả là 0.6. Nói cách khác, cho mỗi lần kiểm định (với dữ liệu null) so sánh AMVR của cả 12 tháng với trung bình, sẽ có 0.6 tháng cho kết quả có ý nghĩa. (rõ ràng không thể có ít hơn 1 tháng, nhưng khi lặp lại nhiều lần, con số sẽ tiến gần giá trị này).
Xuyên suốt công việc này, ta nhấn mạnh tầm quan trọng của việc rõ ràng với câu hỏi cần trả lời. Vấn đề với các kỳ vọng vừa tính là ta đã giả định đang kiểm định kết quả có ý nghĩa so với cả 12 tháng! Vì thế xác suất thấy ít nhất một dương tính giả cao đến 46%.
Đây là ví dụ về sai lệch so sánh nhiều lần khi ta mở rộng không gian tìm kiếm và tăng khả năng tìm được kết quả có ý nghĩa. Điều này nguy hiểm vì có thể bị lạm dụng để chọn lọc các tham số mô hình nhằm đạt p-value ‘mong muốn’.
Đây là tinh túy của p-Hacking
Để minh họa tác động của p-hacking và cách giảm tính đa trị, ta cần hiểu sự khác biệt tinh tế nhưng quan trọng giữa 2 câu hỏi sau:
“Xác suất để tháng 12 trông cực trị như vậy do ngẫu nhiên là bao nhiêu?”
“Xác suất để bất kỳ tháng nào trông cực trị như vậy do ngẫu nhiên là bao nhiêu?”
Vẻ đẹp của mô phỏng nằm ở sự đơn giản. Đoạn mã sau là tất cả những gì cần để tính p-value trả lời câu hỏi thứ nhất. Ta chỉ việc đếm bao nhiêu giá trị trong tập dữ liệu dùng toàn bộ 12000 độ lệch AMVR (12 tháng x 1000 lần thử) lớn hơn giá trị tháng 12 quan sát được. Ta thu được p-value 4.4%, gần ngưỡng 5% có phần tùy tiện của chúng ta, nhưng vẫn có ý nghĩa.
#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 )Được hỗ trợ bởi
Để trả lời câu hỏi thứ hai và tránh tính đa trị, thay vì so sánh kết quả với phân phối tạo từ tất cả 12000 độ lệch AMVR, ta chỉ xét giá trị cao nhất từ mỗi trong số 1000 lần thử AMVR tuyệt đối. Điều này cho p-value 23%, rõ ràng là không có ý nghĩa!
#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 )Được hỗ trợ bởi
Giờ khi đã có kết quả cuối cùng, hãy vẽ các phân phối này để cho thấy tác động của p-Hacking và kết quả phân tích một cách trực quan:
1. Đầu tiên dùng np.quantile() để tìm mức ý nghĩa 5%, ta sẽ vẽ bằng axvline ở các biểu đồ dưới.
2. Tiếp theo định nghĩa 4 ô đồ thị con. Ở đây plt.subplots trả về 'fig' = hình và 'ax1,ax2,ax3,ax4' = biến đại diện cho 4 đồ thị con.
3. #plot 1 hiển thị cột thứ nhất. Ta chỉ việc vẽ biểu đồ tần suất của dữ liệu 'abs_all_months' và chọn type='bar'. Với các biểu đồ dưới, ta đặt cumulative='True', điều này cho CDF thay vì PDF. Bins xác định số cột cho mỗi khoảng giá trị, 30 là con số hợp lý cho dữ liệu của ta. Sau đó định dạng biểu đồ bằng các kỹ thuật đã học, như axvline để vẽ các mức ý nghĩa và giá trị quan sát.
Cột bên trái là dữ liệu trả lời câu hỏi 1 và cột bên phải là câu hỏi 2. Hàng trên là phân phối xác suất, hàng dưới là CDF. Đường gạch đỏ là mức ý nghĩa 5% mà ta đã chọn tùy ý. Đường xanh là giá trị AMVR tháng 12 cực trị gốc là 1.23.
Biểu đồ bên trái cho thấy giá trị tháng 12 gốc có ý nghĩa ở mức 5%, nhưng chỉ vừa đủ! Tuy nhiên, khi tính đến sai lệch so sánh nhiều lần, ở biểu đồ bên phải ngưỡng ý nghĩa tăng từ khoảng 1.2 (abs AMVR) lên khoảng 1.6 (xem đường đỏ).
Khi tính đến sai lệch so sánh nhiều lần, giá trị 1.23 của tháng 12 không còn có ý nghĩa nữa!
Bằng cách cân nhắc kỹ câu hỏi cần trả lời và tránh sai lệch so sánh nhiều lần, ta đã tránh p-hacking mô hình và tránh chỉ ra kết quả có ý nghĩa khi thực tế không có.
Để tìm hiểu thêm về p-hacking và cách bị lạm dụng để kể một câu chuyện nhất định về dữ liệu, hãy xem ứng dụng tương tác tuyệt vời từ FiveThirtyEight
Kết luận
Chúng ta đã học rằng kiểm định giả thuyết không đáng sợ như tưởng tượng. Chỉ cần làm theo 3 bước trên để xây dựng mô hình cho bất kỳ loại dữ liệu hay thống kê kiểm định nào.
Chúng ta đã cho thấy việc đặt đúng câu hỏi là tối quan trọng cho phân tích khoa học. Một thay đổi nhỏ trong cách diễn đạt có thể dẫn đến mô hình rất khác và kết quả rất khác.
Hy vọng sức mạnh của các chức năng Python trung cấp và khả năng hỗ trợ bạn thực hiện kiểm định thống kê đã trở nên rõ ràng! Chỉ với vài dòng mã, chúng ta đã xây dựng và kiểm định một hiện tượng thực tế và rút ra kết luận có thể hành động.
Chúng ta đã thảo luận tầm quan trọng của việc nhận biết và hiệu chỉnh sai lệch so sánh nhiều lần và tránh cạm bẫy p-hacking, đồng thời cho thấy một kết quả có vẻ có ý nghĩa có thể trở thành không có ý nghĩa.
Với ngày càng nhiều ‘dữ liệu lớn’ cùng áp lực học thuật phải tạo ra nghiên cứu với phát hiện ‘mới lạ’ hoặc áp lực chính trị để cho thấy một kết quả là ‘có ý nghĩa’, cám dỗ p-hacking ngày càng tăng. Bằng cách học cách nhận ra khi nào ta mắc phải và hiệu chỉnh tương ứng, chúng ta có thể trở thành những nhà nghiên cứu tốt hơn và cuối cùng tạo ra kết quả khoa học chính xác hơn, từ đó có thể hành động!
Ghi chú của tác giả: Kết quả của chúng ta hơi khác với nghiên cứu gốc của Winton, một phần do bộ dữ liệu khác (32 năm so với 87 năm) và họ chọn tháng 10 là tháng đáng chú ý trong khi chúng ta là tháng 12. Ngoài ra, họ dùng phương pháp ‘dữ liệu mô phỏng’ không công khai, trong khi chúng ta đã nêu rõ qua mã phương pháp tạo dữ liệu đó. Chúng ta đã đưa ra một số giả định mô hình xuyên suốt, những điều này đã được nêu rõ và có thể thấy trong mã. Những lựa chọn thiết kế và mô hình hóa là một phần của quy trình khoa học, miễn là được trình bày minh bạch, phân tích vẫn có giá trị.
Rèn kỹ năng Python để thăng tiến trong ngành tài chính. Học cách làm việc với danh sách, mảng và trực quan hóa dữ liệu để thành thạo phân tích tài chính.
%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)
%22%3E%0A%20%20%20%20%20%20%20%20%3Crect%20fill%3D%22%23EFEBE4%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%23DD3400%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%23FFEC3C%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%23FF931E%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%23DD3400%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%23FFEC3C%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%23FF931E%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)
Matt Crabtree · 10 min
"Trong 32 năm qua, tháng 10 là tháng biến động nhất tính trung bình đối với S&P500 và tháng 12 là tháng ít biến động nhất".
Trong hướng dẫn này, chúng ta sẽ dùng Python để đi qua toàn bộ quá trình phân tích và kiểm định hiện tượng này nhằm xác định liệu nó có ý nghĩa thống kê hay không.
Chúng ta sẽ sử dụng các chức năng chuỗi thời gian của pandas, như Resample, để biến đổi dữ liệu giá cổ phiếu thô thành định dạng làm dữ liệu trở nên trực quan.
Sau đó, chúng ta sẽ tìm hiểu cách thực hiện kiểm định giả thuyết bằng mô phỏng trong Python thay vì dùng các công thức rườm rà.
Cuối cùng, chúng ta sẽ bàn về một vấn đề lớn trong phân tích thống kê, hiện tượng sai lệch so sánh nhiều lần (Multiple Comparison Bias), cách xử lý và minh họa trực quan qua matplotlib tác động của 'p-Hacking'.
Mục tiêu của chúng ta:
Dữ liệu của chúng ta:
Chúng ta sẽ sử dụng dữ liệu S&P500 theo ngày cho phân tích này, cụ thể là giá đóng cửa thô hằng ngày từ 1986 đến 2018 (khá khó tìm, nên tôi đã công khai).
Cảm hứng cho bài viết đến từ Winton, chúng ta sẽ tái hiện lại ở đây, dù chỉ với 32 năm dữ liệu so với 87 năm của họ.
Xử lý dữ liệu với Pandas
Để trả lời câu hỏi liệu biến động cực đoan ở một số tháng có thực sự đáng kể và do đó có khả năng tiếp diễn hay không, chúng ta cần biến đổi 32 năm dữ liệu giá thành một định dạng thể hiện rõ hiện tượng đang khảo sát.
Định dạng lựa chọn là bảng xếp hạng biến động trung bình theo tháng (AMVR).
Đoạn mã sau cho thấy cách đưa dữ liệu giá thô về định dạng này. Bắt đầu thôi!
Đầu tiên là các import tiêu chuẩn. (matplotlib.patches cho phép chúng ta kiểm soát việc tạo kiểu cho từng cột trong biểu đồ cột)
Một mẹo hay để biểu đồ hiển thị toàn chiều rộng trong Jupyter notebooks:
Cái này tùy chọn - hãy thử và xem điều gì xảy ra!
Import dữ liệu bằng phương thức read_csv nhận một đường dẫn, ở đây là url đến dữ liệu. Ta cũng chỉ định dùng 'date' làm chỉ mục và tự động phân tích cú pháp ngày tháng từ văn bản cho phù hợp nhất có thể.
Điều này cho ta giá đóng cửa thô chưa điều chỉnh của S&P500 (SPX). Giờ ta cần chuyển giá thô sang lợi nhuận % theo ngày. Có hai lựa chọn:
Với mục đích của chúng ta, sẽ dùng cách 2, vì pandas có thể xử lý tức thì tập dữ liệu cỡ này (hơn 8000 giá trị). Ta cũng sẽ làm sạch dữ liệu bằng cách bỏ giá trị đầu tiên trở thành 'NaN' do không có thay đổi giá so với ngày trước đó. pct_change() nhận tham số tùy chọn 'periods' để đổi khoảng dịch; ta giữ mặc định là 1.
Bước tiếp theo là lấy các giá trị % thay đổi theo ngày và chuyển chúng thành 'biến động hàng tháng đã quy đổi theo năm'. Dòng đầu tiên dưới đây cho thấy ta có thể thực hiện chuyển đổi này chỉ bằng một dòng mã. Đó là sức mạnh của Pandas! Nhưng hãy bóc tách từng bước.
1. Để có mnthly_annu, đầu tiên ta dùng phương thức 'resample' trên lợi nhuận ngày. Resample cho phép thay đổi tần suất kỳ dữ liệu. Nó nhận một 'chuỗi ký hiệu tần suất', nói đơn giản là một chữ cái tương ứng tần suất mong muốn. Bao gồm:
Vì muốn lấy mẫu lại từ lợi nhuận ngày sang tháng, ta truyền 'M' làm tham số.
2. Tiếp theo cần quyết định cách tính giá trị hàng tháng mới, ví dụ có thể cộng dồn, nhân, v.v. Với phân tích của ta, cần thước đo biến động và độ lệch chuẩn phù hợp, nên ta nối thêm std() vào dữ liệu đã resample để lấy biến động tháng.
3. Bước cuối là quy đổi theo năm con số này. Ta làm đơn giản bằng cách nhân với căn bậc hai của 12, với 12 là số kỳ (tháng) trong năm. Điều này cho ta các giá trị biến động tháng quy đổi theo năm cần thiết.
Thực hành tốt trong phân tích dữ liệu là trực quan hóa dữ liệu khi tiến hành phân tích. Hãy xem dữ liệu biến động đã quy đổi theo năm của chúng ta.
Biểu đồ dưới cho thấy rõ các sự kiện lớn trên thị trường, như Black Monday và Khủng hoảng Tài chính 2008. Phương thức 'axvspan' của matplotlib cho phép thêm các vạch dọc (axhspan là phương thức tương ứng cho phương ngang), nhận 'xmin' và 'xmax' làm tham số, xác định bề rộng của hình chữ nhật tạo ra; vì trục x là chỉ mục datetime, ta có thể truyền trực tiếp các năm muốn làm nổi bật. Tham số alpha điều chỉnh độ trong suốt của khối màu để vẫn thấy biểu đồ bên dưới.
Thuộc tính mpatches cho phép tạo một hộp chú giải tùy chỉnh. Ta làm bằng cách định nghĩa biến 'labs' với màu, alpha và văn bản mong muốn, sau đó truyền biến này vào tham số 'handles' trong plt.legend để hiển thị chú giải.
Vậy là chúng ta đã thấy một trong những phương thức mạnh mẽ của pandas, resample. Giờ hãy dùng một phương thức khác, 'groupby'. Chúng ta cần đi từ các giá trị biến động tháng đã quy đổi theo năm sang thước đo AMVR mong muốn. Lại một lần nữa, ta có thể đạt được chỉ với vài dòng mã.
1. Đầu tiên, áp dụng 'groupby' lên mnthly_annu. Groupby cần một tham số xác định cách nhóm dữ liệu. Có thể là một hàm hoặc như ở đây, một Series. Ta truyền 'mnthly_annu.index.year', đơn giản là thuộc tính năm trên chỉ mục datetime của mnthly_annu. Điều này sẽ nhóm các giá trị biến động tháng cho từng năm trong bộ dữ liệu 32 năm.
2. Tiếp theo, áp dụng phương thức 'rank', sắp xếp dữ liệu theo thứ tự tăng dần.
3. Cuối cùng, lặp lại và lấy trung bình theo tất cả các năm cho từng tháng. Điều này cho ra các giá trị AMVR cuối cùng!
Điều này cho ta Bảng xếp hạng Biến động Trung bình theo Tháng cuối cùng. Về số liệu, ta thấy tháng 10 là cao nhất và tháng 12 là thấp nhất.
Chọn trực quan hóa phù hợp cho dữ liệu là quan trọng vì chúng ta muốn kể câu chuyện dữ liệu qua những biểu đồ tác động nhất. Với phân tích này, ta muốn thể hiện rõ mức biến động cao nhất và thấp nhất. Ta sẽ dùng biểu đồ cột của matplotlib và chú thích, tô màu phù hợp để tăng tác động.
1. Bằng cách lập chỉ mục b_plot (b_plot[9]) với cột mong muốn, ta có thể đặt màu để làm nổi bật giá trị cao nhất và thấp nhất.
2. Để thêm giá trị AMVR lên cột tương ứng, ta duyệt liệt kê qua từng giá trị AMVR (biến 'final') và làm tròn đến 2 chữ số thập phân. 'i' và 'v' lần lượt là 'index' và 'value' của biến 'final', nên 'i' sẽ từ 1-12, và 'value' là giá trị AMVR làm tròn 2 chữ số. Ta dùng các biến này trong plt.text(). Tham số thứ nhất của plt.text() là vị trí trục x cho nhãn, nên ta truyền 'i' cho mỗi cột và dịch 0.8 để căn giữa. Tham số thứ hai là chuỗi biểu diễn giá trị 'final'.
3. Để hiển thị giá trị trung bình, ta dùng axhline và dùng 'label' để thêm chú giải ở góc trên bên phải. Chỉ cần truyền giá trị trục y, tức giá trị trung bình, rồi thêm định dạng.
Biểu đồ này là trực quan tốt nhất cho hiện tượng đang khảo sát. Ta thấy rõ tháng 10 biến động mạnh nhất và tháng 12 thấp nhất. Quan trọng là ta cũng thấy tháng 12 là giá trị 'cực trị' hơn, theo giá trị tuyệt đối. Điều này sẽ quan trọng cho phân tích ở phần tiếp theo.
Vậy là xong phần dữ liệu, giờ đến Kiểm định giả thuyết…
Kiểm định giả thuyết: Câu hỏi là gì?
Kiểm định giả thuyết là một trong những kỹ thuật nền tảng của khoa học dữ liệu, nhưng cũng là một trong những thứ gây e ngại và dễ bị hiểu sai. Căn nguyên nằm ở cách dạy trong môn Thống kê cơ bản, nơi ta được bảo:
"thực hiện t-test, đó là kiểm định một phía hay hai phía? Chọn thống kê kiểm định phù hợp như t-test của Welch, tính bậc tự do, tính điểm t, tra giá trị tới hạn trong bảng, so sánh giá trị tới hạn với thống kê t ……"
Dễ hiểu khi điều này dẫn đến bối rối về việc nên dùng kiểm định nào và tiến hành ra sao. Tuy nhiên, tất cả các kỹ thuật thống kê cổ điển để thực hiện kiểm định giả thuyết được phát triển vào thời điểm khả năng tính toán rất hạn chế và chỉ là các nghiệm đóng để tính p-value, hết! Nhưng kèm thêm sự phức tạp là phải chọn đúng công thức cho tình huống cụ thể, do các giả định hạn chế và đôi khi mù mờ.
Nhưng hãy mừng lên!
Có cách tốt hơn. Mô phỏng.
Để hiểu mô phỏng giúp gì, hãy nhắc lại kiểm định giả thuyết là gì:
Chúng ta muốn kiểm tra “liệu hiệu ứng quan sát được trong dữ liệu là thật hay có thể xảy ra đơn thuần do ngẫu nhiên” và để thực hiện, ta làm như sau:
Chỉ có một bài kiểm định! Bài blog hay của Allen Downey có thêm chi tiết về kiểm định giả thuyết
Sức mạnh thực sự của mô phỏng là buộc chúng ta phải diễn đạt rõ ràng các giả định mô hình qua mã. Trong khi kỹ thuật cổ điển có thể là một ‘hộp đen’ về các giả định của chúng.
Ví dụ dưới đây: Biểu đồ trái là dữ liệu thật và hiệu ứng quan sát với một xác suất nhất định (màu xanh lá). Biểu đồ phải là dữ liệu null mô phỏng với ghi nhận khi nào hiệu ứng quan sát xuất hiện do ngẫu nhiên (màu đỏ). Đây là nền tảng của kiểm định giả thuyết: xác suất thấy hiệu ứng quan sát trên dữ liệu null là bao nhiêu.
Phần quan trọng nhất của kiểm định giả thuyết là rõ ràng về câu hỏi cần trả lời. Ở đây chúng ta hỏi:
“Giá trị cực trị nhất có thể xảy ra do ngẫu nhiên không?”
Giá trị cực trị nhất được định nghĩa là độ lệch tuyệt đối của AMVR so với trung bình lớn nhất. Câu hỏi này tạo nên giả thuyết không của chúng ta.
Trong dữ liệu, giá trị cực trị nhất là giá trị tháng 12 (1.23) chứ không phải tháng 10 (1.08), vì ta xét độ lệch tuyệt đối lớn nhất so với trung bình, không chỉ là biến động cao nhất.
1. Để có độ lệch tuyệt đối, chỉ cần lấy mỗi giá trị trừ đi trung bình và áp dụng hàm abs().
2. Dùng sort_values() để sắp xếp kết quả từ nhỏ đến lớn. Sau đó chọn ra 2 giá trị lớn nhất, Tháng 10 & Tháng 12.
Mô phỏng
Giờ ta đã biết câu hỏi cần hỏi, cần xây dựng ‘Mô hình Null’.
Có một số lựa chọn:
Lấy mẫu là một chủ đề lớn, nhưng ta sẽ tập trung vào một kỹ thuật, hoán vị hay xáo trộn.
Để có mô hình null mong muốn, ta cần dựng một tập dữ liệu không có tính mùa vụ. Nếu giả thuyết không đúng, rằng không có mùa vụ trong dữ liệu và hiệu ứng quan sát chỉ do ngẫu nhiên, thì nhãn cho từng tháng (Jan, Feb, v.v.) là vô nghĩa và do đó ta có thể xáo trộn dữ liệu lặp lại để xây dựng thứ mà thống kê cổ điển gọi là ‘phân phối lấy mẫu của thống kê kiểm định dưới giả thuyết không’. Điều này khiến hiện tượng quan sát (giá trị tháng 12 cực trị) có khả năng như nhau cho mọi tháng, đúng như yêu cầu của mô hình null.
Để chứng minh sức mạnh của kỹ thuật mô phỏng với năng lực tính toán hiện đại, mã trong ví dụ này thực sự hoán vị dữ liệu giá theo ngày, đòi hỏi nhiều sức mạnh xử lý hơn, nhưng vẫn hoàn tất trong vài giây trên CPU hiện đại.
Lưu ý: xáo trộn nhãn theo ngày hoặc theo tháng đều cho tập dữ liệu null mong muốn trong trường hợp này.
Một tài nguyên tuyệt vời để học về lấy mẫu đến từ Julian Simon.
Lưu ý: Cách ta xây dựng kiểm định tương đương kiểm định hai phía theo phương pháp ‘cổ điển’, như t-test của Welch hay ANOVA, v.v., vì ta quan tâm đến giá trị cực trị nhất, ở trên hay dưới trung bình.
Những quyết định này là lựa chọn thiết kế, và ta có quyền tự do đó vì Mô hình Null chỉ là mô hình! Điều này nghĩa là ta có thể chỉ định tham số theo ý, điều mấu chốt là phải thật rõ ràng câu hỏi cần trả lời.
Vậy ta muốn dùng Python để mô phỏng nhiều dữ liệu nhằm tạo tập dữ liệu null. Để làm điều này, ta sẽ mô phỏng 1000 bộ gồm 12 AMVR, hoán vị nhãn 'date' mỗi lần để xây dựng phân phối lấy mẫu. Kết quả từ mã này được đưa ở phần p-hacking bên dưới.
1. Đầu tiên ta định nghĩa vài cấu trúc chứa kết quả: sẽ dùng một pd.Dataframe() và một mảng đơn [ ].
2. Tiếp theo, định nghĩa một biến đếm và đặt 0, rồi bắt đầu vòng lặp for 1000 lần để tạo dữ liệu mô phỏng.
3. Dòng đầu trong vòng lặp lấy lợi nhuận ngày gốc (từ đầu hướng dẫn) và dùng phương thức sample() của pandas. Phương thức này lấy mẫu ngẫu nhiên từ dữ liệu lợi nhuận ngày một số lần nhất định. Ở đây, ta có 8191 điểm dữ liệu (252 ngày giao dịch trong 32 năm). Ta cũng bỏ chỉ mục bằng reset_index(), để có thể thêm chỉ mục mới. Cần làm vậy vì chỉ mục ngày gốc sẽ bị xáo trộn cùng dữ liệu, trong khi ta muốn 'xáo trộn' chỉ dữ liệu.
4. Dòng tiếp theo thêm chỉ mục ngày mới bằng cách gán lại chỉ mục của daily_ret_shuffle với một pd.bdate_range có độ dài bằng dữ liệu gốc. Ta dùng bdate_range thay vì date_range vì cần ngày làm việc (5, Thứ Hai - Thứ Sáu) thay vì 7 ngày một tuần.
5. Giờ đã tạo dữ liệu 'xáo trộn', ta thực hiện lại quá trình xử lý như lúc đầu để tạo giá trị AMVR. Ba dòng mã tiếp theo làm việc này.
6. Khi đã có giá trị AMVR mô phỏng, ta cần nối chúng vào dataframe lưu trữ toàn bộ 1000 lượt mô phỏng. pd.concat nhận dataframe hiện tại và nối từng dataframe tiếp theo vào cuối. Ta chọn axis 1, để có các cột dữ liệu.
7. Maximonth là biến chỉ lưu giá trị cao nhất từ mỗi lần mô phỏng (sẽ dùng sau trong phân tích).
8. Vì phân tích cần giá trị AMVR tuyệt đối, 3 dòng mã tiếp theo lấy toàn bộ 1000 lượt AMVR và làm phẳng thành một mảng lớn. Sau đó lấy trung bình và trừ từng giá trị cho trung bình này rồi lấy giá trị tuyệt đối.
9. Ta cũng làm tương tự cho dữ liệu 'chỉ cao nhất'. Lưu ý: ở đây dùng list thay vì dataframe nên có thể dùng list comprehension để tính abs(AMVR) cho từng giá trị cao nhất.
Giờ ta đã có mọi dữ liệu cần để hoàn tất phân tích, cả quan sát gốc lẫn tập dữ liệu mô phỏng mới.
p-Hacking
Đây là phần thú vị. Ta đã xây dựng giả thuyết để kiểm định, đã tạo dữ liệu mô phỏng bằng cách xáo trộn nhãn ‘date’, giờ cần thực hiện kiểm định để tìm xác suất quan sát được kết quả đáng kể như kết quả tháng 12 với giả thuyết không (không có mùa vụ) là đúng.
Trước khi kiểm định, hãy đặt kỳ vọng.
Xác suất thấy ít nhất một kết quả có ý nghĩa với mức ý nghĩa 5% là bao nhiêu?
= 1-p(không có ý nghĩa)
= 1-(1–0.05)¹²
= 0.46
vậy có 46% khả năng thấy ít nhất một tháng có kết quả có ý nghĩa, giả thuyết không là đúng.
Giờ hãy hỏi, với từng kiểm định riêng lẻ (so sánh độ lệch tuyệt đối AMVR của 12 tháng so với trung bình) chúng ta kỳ vọng thấy bao nhiêu giá trị có ý nghĩa trong dữ liệu ngẫu nhiên, không mùa vụ?
12 x 0.05 = 0.6
Vậy với mức ý nghĩa 0.05, ta kỳ vọng tỷ lệ dương tính giả là 0.6. Nói cách khác, cho mỗi lần kiểm định (với dữ liệu null) so sánh AMVR của cả 12 tháng với trung bình, sẽ có 0.6 tháng cho kết quả có ý nghĩa. (rõ ràng không thể có ít hơn 1 tháng, nhưng khi lặp lại nhiều lần, con số sẽ tiến gần giá trị này).
Xuyên suốt công việc này, ta nhấn mạnh tầm quan trọng của việc rõ ràng với câu hỏi cần trả lời. Vấn đề với các kỳ vọng vừa tính là ta đã giả định đang kiểm định kết quả có ý nghĩa so với cả 12 tháng! Vì thế xác suất thấy ít nhất một dương tính giả cao đến 46%.
Đây là ví dụ về sai lệch so sánh nhiều lần khi ta mở rộng không gian tìm kiếm và tăng khả năng tìm được kết quả có ý nghĩa. Điều này nguy hiểm vì có thể bị lạm dụng để chọn lọc các tham số mô hình nhằm đạt p-value ‘mong muốn’.
Để minh họa tác động của p-hacking và cách giảm tính đa trị, ta cần hiểu sự khác biệt tinh tế nhưng quan trọng giữa 2 câu hỏi sau:
Vẻ đẹp của mô phỏng nằm ở sự đơn giản. Đoạn mã sau là tất cả những gì cần để tính p-value trả lời câu hỏi thứ nhất. Ta chỉ việc đếm bao nhiêu giá trị trong tập dữ liệu dùng toàn bộ 12000 độ lệch AMVR (12 tháng x 1000 lần thử) lớn hơn giá trị tháng 12 quan sát được. Ta thu được p-value 4.4%, gần ngưỡng 5% có phần tùy tiện của chúng ta, nhưng vẫn có ý nghĩa.
Để trả lời câu hỏi thứ hai và tránh tính đa trị, thay vì so sánh kết quả với phân phối tạo từ tất cả 12000 độ lệch AMVR, ta chỉ xét giá trị cao nhất từ mỗi trong số 1000 lần thử AMVR tuyệt đối. Điều này cho p-value 23%, rõ ràng là không có ý nghĩa!
Giờ khi đã có kết quả cuối cùng, hãy vẽ các phân phối này để cho thấy tác động của p-Hacking và kết quả phân tích một cách trực quan:
1. Đầu tiên dùng np.quantile() để tìm mức ý nghĩa 5%, ta sẽ vẽ bằng axvline ở các biểu đồ dưới.
2. Tiếp theo định nghĩa 4 ô đồ thị con. Ở đây plt.subplots trả về 'fig' = hình và 'ax1,ax2,ax3,ax4' = biến đại diện cho 4 đồ thị con.
3. #plot 1 hiển thị cột thứ nhất. Ta chỉ việc vẽ biểu đồ tần suất của dữ liệu 'abs_all_months' và chọn type='bar'. Với các biểu đồ dưới, ta đặt cumulative='True', điều này cho CDF thay vì PDF. Bins xác định số cột cho mỗi khoảng giá trị, 30 là con số hợp lý cho dữ liệu của ta. Sau đó định dạng biểu đồ bằng các kỹ thuật đã học, như axvline để vẽ các mức ý nghĩa và giá trị quan sát.
Cột bên trái là dữ liệu trả lời câu hỏi 1 và cột bên phải là câu hỏi 2. Hàng trên là phân phối xác suất, hàng dưới là CDF. Đường gạch đỏ là mức ý nghĩa 5% mà ta đã chọn tùy ý. Đường xanh là giá trị AMVR tháng 12 cực trị gốc là 1.23.
Biểu đồ bên trái cho thấy giá trị tháng 12 gốc có ý nghĩa ở mức 5%, nhưng chỉ vừa đủ! Tuy nhiên, khi tính đến sai lệch so sánh nhiều lần, ở biểu đồ bên phải ngưỡng ý nghĩa tăng từ khoảng 1.2 (abs AMVR) lên khoảng 1.6 (xem đường đỏ).
Bằng cách cân nhắc kỹ câu hỏi cần trả lời và tránh sai lệch so sánh nhiều lần, ta đã tránh p-hacking mô hình và tránh chỉ ra kết quả có ý nghĩa khi thực tế không có.
Để tìm hiểu thêm về p-hacking và cách bị lạm dụng để kể một câu chuyện nhất định về dữ liệu, hãy xem ứng dụng tương tác tuyệt vời từ FiveThirtyEight
Kết luận
Ghi chú của tác giả: Kết quả của chúng ta hơi khác với nghiên cứu gốc của Winton, một phần do bộ dữ liệu khác (32 năm so với 87 năm) và họ chọn tháng 10 là tháng đáng chú ý trong khi chúng ta là tháng 12. Ngoài ra, họ dùng phương pháp ‘dữ liệu mô phỏng’ không công khai, trong khi chúng ta đã nêu rõ qua mã phương pháp tạo dữ liệu đó. Chúng ta đã đưa ra một số giả định mô hình xuyên suốt, những điều này đã được nêu rõ và có thể thấy trong mã. Những lựa chọn thiết kế và mô hình hóa là một phần của quy trình khoa học, miễn là được trình bày minh bạch, phân tích vẫn có giá trị.
Nếu bạn muốn học thêm về tài chính trong Python, hãy học khóa Intro to Python for Finance của DataCamp và xem Python For Finance Tutorial: Algorithmic Trading.
Theo dõi tôi trên twitter.com/pdquant để biết thêm!