-
Notifications
You must be signed in to change notification settings - Fork 0
/
main.py
129 lines (92 loc) · 4.45 KB
/
main.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
import pickle
from scipy.stats import t, norm
import matplotlib.pyplot as plt
import numpy as np
from eventstudystatistics import adjBMP, adjBMP_daily, grank
def calculate_coefficients(X,Y):
# implement ordinary least squares in numpy
# add a constant to the X matrix
X = np.c_[np.ones(X.shape[0]), X]
# calculate the coefficients
beta = np.linalg.inv(X.T @ X) @ X.T @ Y
# calculate the residuals
eps = Y - X @ beta
return beta[0], beta[1], eps
# load pickle file from tests
with open("tests/params_adjbmp_grank.pkl", "rb") as f:
params = pickle.load(f)
test1 = grank(*params)
test2 = adjBMP(*params)
CAR_period = [0, 40] # including both edges, CAREFUL, this is not like python indexing, this is including the right side.
n_events = 120
length_event_window = 41 # L2
length_estimation_window = 100 # L1
event_day = 20 # within the event window index 20 is the event day
# single example for all tests:
event_window_market_return = np.random.normal(0, 0.1, (n_events, length_event_window))
event_window_company_return = np.random.normal(0, 0.05, (n_events, length_event_window)) + event_window_market_return
estimation_window_market_return = np.random.normal(0, 0.1, (n_events, length_estimation_window))
estimation_window_company_return = np.random.normal(0, 0.05, (n_events, length_estimation_window)) + estimation_window_market_return
AR_ = []
eps_ = []
print("calculate abnormal returns...")
for i in range(n_events):
alpha, beta, eps = calculate_coefficients(estimation_window_market_return[i, :],
estimation_window_company_return[i, :])
## Calculate the abnormal returns
abnormal_return = event_window_company_return[i, :] - alpha - beta * event_window_market_return[i, :]
AR_.append(abnormal_return)
eps_.append(eps)
print("Done calculating abnormal returns")
AR = np.asarray(AR_)
eps = np.asarray(eps_)
test_res = adjBMP_daily(AR, eps, estimation_window_market_return, event_window_market_return, event_day)
print(test_res)
test_res = adjBMP(AR, eps, estimation_window_market_return, event_window_market_return, CAR_period)
print(test_res)
test_res2 = grank(AR, eps, estimation_window_market_return, event_window_market_return, CAR_period)
print(test_res2)
### SIMULATION:
### test if there are no missing values in the data, expecting high p value
adjBMP_results = []
grank_results = []
np.random.seed(3)
J = 5000
for j in range(J):
print(j)
event_window_market_return = np.random.normal(0, 0.1, (n_events, length_event_window))
event_window_company_return = np.random.normal(0, 0.05, (n_events, length_event_window)) + event_window_market_return
estimation_window_market_return = np.random.normal(0, 0.1, (n_events, length_estimation_window))
estimation_window_company_return = np.random.normal(0, 0.05, (n_events, length_estimation_window)) + estimation_window_market_return
AR_ = []
eps_ = []
print("calculate AR...")
for i in range(n_events):
alpha, beta, eps = calculate_coefficients(estimation_window_market_return[i, :],
estimation_window_company_return[i, :])
## Calculate the abnormal returns
abnormal_return = event_window_company_return[i, :] - alpha - beta * event_window_market_return[i, :]
AR_.append(abnormal_return)
eps_.append(eps)
print("Done calculating AR")
AR = np.asarray(AR_)
eps = np.asarray(eps_)
#test_res = adjBMP(AR, eps, estimation_window_market_return, event_window_market_return, CAR_period, adjustment=False)
#adjBMP_results.append(test_res)
#test_res_daily = adjBMP_daily(AR, eps, estimation_window_market_return, event_window_market_return, event_day, adjustment=False)
#adjBMP_results.append(test_res_daily)
test_res2 = grank(AR, eps, estimation_window_market_return, event_window_market_return, CAR_period, True)
grank_results.append(test_res2)
adj_BMP_z_stat = np.asarray([res.statistic for res in adjBMP_results])
grank_t_stat = np.asarray([res.statistic for res in grank_results])
# histogram of the statistics
plt.hist(adj_BMP_z_stat, bins=200, density=True)
# add a standard normal distribution
x = np.linspace(-5, 5, length_estimation_window)
plt.plot(x, norm.pdf(x))
plt.show()
plt.hist(grank_t_stat, bins=200, density=True)
# add a student t distribution of 99 degrees of freedom
x = np.linspace(-5, 5, length_estimation_window)
plt.plot(x, t.pdf(x, length_estimation_window-1))
plt.show()