<jupyter_start><jupyter_text>Load previous sesssion<jupyter_code>import dill
dill.load_session('DJIA_univHP_env_September_version2.db')<jupyter_output><empty_output><jupyter_text>Load the data discussed and analysed in https://methods.sagepub.com/dataset/howtoguide/multivariate-hawkes-in-djia-2018i100<jupyter_code># Read data
# Import the data as a dataframe (2D data structure with labelled axes)

df = pd.read_csv('financial_data/dataset-djia-2018-subset2.csv')
dates = pd.to_datetime(df['Date']) # set dates as the Date of closing price column


del df['Date'] # delete Date column

# Fill missing values
df.ffill(inplace=True)
df
# Google and amazon are not in DJIA list

ten_companies = ['AAPL','MSFT', 'JPM', 'GS', 'PFE', 'MRK','NKE', 'HD',  'GOOGL','AMZN']


ticker_id = [] # collect integer marker of each stock, index correspodns to stock AABA, 1 correponds to stock AAPL
ticker = []
ticker_dict = {} # collect key-value pairs where key represents tikcer of a stock and it's corresponding value an integer marker


for u,col in enumerate(df[ten_companies]):
    ticker_dict[col] = u
    ticker_id += [u]
    ticker += [col]<jupyter_output><empty_output><jupyter_text>Collect relevant data such as event times (largest $10 \%$ price drop, only include days were return was below 10% quantile), markers,...<jupyter_code># Collect event times

t_i=[]
u_i=[] # collect marker for each event times

ticker_id = [] # collect integer marker of each stock, index correspodns to stock AABA, 1 correponds to stock AAPL
ticker = []
ticker_dict = {} # collect key-value pairs where key represents tikcer of a stock and it's corresponding value an integer marker


for u,col in enumerate(df[ten_companies]):
    Tdiff = df[col].diff()
    timestamps=dates[Tdiff<=Tdiff.quantile(0.1)] # return timestamps where subsequent price difference was less than 10% quantile
    t_i.extend((timestamps - pd.Timestamp(2006,1,3)).dt.days.astype(float)) # measure time 
    u_i.extend(np.repeat(u,len(timestamps)))
    ticker_dict[col] = u
    ticker_id += [u]
    ticker += [col]

t_i = np.array(t_i)
#t_i=np.array(t_i)+np.random.rand(len(t_i))
u_i=np.array(u_i)
perm = np.argsort(t_i)
t_i = t_i[perm] # superposed timestamps
u_i = u_i[perm]
tstamps_ten_univ = []
tstamps_ten_univ_dates = []

count = 0
for i in ticker:
    tdiff_i = df[i].diff()
    tstamps_ten_univ_dates += [dates[tdiff_i<=tdiff_i.quantile(0.1)]]
    tstamps_i = (tstamps_ten_univ_dates[count] - pd.Timestamp(2006,1,3)).dt.days.astype(float) # change time units from days to months
    tstamps_ten_univ += [tstamps_i]
    count += 1
ts_i_with_dates = []
for i in range(len(ten_companies)):
    ts_i_with_dates += [pd.DataFrame({'Date': tstamps_ten_univ_dates[i].values, 'event time': tstamps_ten_univ[i].values})]
    #ts_i_with_dates[i].to_csv('ten_stocks_event_times/'+ten_companies[i]+'_ts_with_dates.csv', index=None)
avg1 = 0
avg2 = 0

for i in ticker:
    total_events = t_i[np.where(u_i == ticker_dict[i])].shape[0]
    print(f'number of events of stock '+i+':', total_events)
    print(f'average time between each event of \'type\' '+i+':', t_i[np.where(u_i == ticker_dict[i])][-1]/total_events)
    avg1 += total_events
    avg2 += t_i[np.where(u_i == ticker_dict[i])][-1]/total_events
print('')
print(f'Total number of total events for 10 companies: ',avg1)

print(f'Average number of total events for 10 companies: ',avg1/10)
print(f'Average time between each event for 10 companies: ',avg2/10)<jupyter_output>number of events of stock AAPL: 303
average time between each event of 'type' AAPL: 14.448844884488448
number of events of stock MSFT: 305
average time between each event of 'type' MSFT: 14.321311475409836
number of events of stock JPM: 302
average time between each event of 'type' JPM: 14.443708609271523
number of events of stock GS: 303
average time between each event of 'type' GS: 14.415841584158416
number of events of stock PFE: 303
average time between each event of 'type' PFE: 14.41914191419142
number of events of stock MRK: 304
average time between each event of 'type' MRK: 14.351973684210526
number of events of stock NKE: 302
average time between each event of 'type' NKE: 14.490066225165563
number of events of stock HD: 302
average time between each event of 'type' HD: 14.437086092715232
number of events of stock GOOGL: 302
average time between each event of 'type' GOOGL: 14.413907284768213
number of events of stock AMZN: 302
average time between each event of 'type' AMZN: 14.4[...]<jupyter_text>Reshape timestamps and marker data for largest $10 \%$ price-jumps into the form of Steve Morse's EM code, and convert unit of time from days to months, i.e. scale each event times by 30<jupyter_code>data = []
counter = 0
for t,u in zip(t_i, u_i):
    data.append([])
    data[counter] += [t/30, np.array([u], dtype='int')]
    counter += 1


data = np.array(data)
data<jupyter_output><empty_output><jupyter_text>Collect relevant data such as event times (largest 15% price drop, only include days were return was below 15% quantile), markers,...<jupyter_code># Collect event times

t_i=[]
u_i=[] # collect marker for each event times

ticker_id = [] # collect integer marker of each stock, index correspodns to stock AABA, 1 correponds to stock AAPL
ticker = []
ticker_dict = {} # collect key-value pairs where key represents tikcer of a stock and it's corresponding value an integer marker


for u,col in enumerate(df[ten_companies]):
    Tdiff = df[col].diff()
    timestamps=dates[Tdiff<=Tdiff.quantile(0.15)] # return timestamps where subsequent price difference was less than 10% quantile
    t_i.extend((timestamps - pd.Timestamp(2006,1,3)).dt.days.astype(float)) # measure time 
    u_i.extend(np.repeat(u,len(timestamps)))
    ticker_dict[col] = u
    ticker_id += [u]
    ticker += [col]

t_i = np.array(t_i)
#t_i=np.array(t_i)+np.random.rand(len(t_i))
u_i=np.array(u_i)
perm = np.argsort(t_i)
t_i = t_i[perm] # superposed timestamps
u_i = u_i[perm]
tstamps_ten_univ = []
tstamps_ten_univ_dates = []

count = 0
for i in ticker:
    tdiff_i = df[i].diff()
    tstamps_ten_univ_dates += [dates[tdiff_i<=tdiff_i.quantile(0.15)]]
    tstamps_i = (tstamps_ten_univ_dates[count] - pd.Timestamp(2006,1,3)).dt.days.astype(float) # change time units from days to months
    tstamps_ten_univ += [tstamps_i]
    count += 1
ts_i_with_dates = []
for i in range(len(ten_companies)):
    ts_i_with_dates += [pd.DataFrame({'Date': tstamps_ten_univ_dates[i].values, 'event time': tstamps_ten_univ[i].values})]
    #ts_i_with_dates[i].to_csv('ten_stocks_event_times/'+ten_companies[i]+'_ts_with_dates.csv', index=None)
print(ts_i_with_dates[0]['event time'])
avg1 = 0
avg2 = 0

for i in ticker:
    total_events = t_i[np.where(u_i == ticker_dict[i])].shape[0]
    print(f'number of events of stock '+i+':', total_events)
    print(f'average time between each event of \'type\' '+i+':', t_i[np.where(u_i == ticker_dict[i])][-1]/total_events)
    avg1 += total_events
    avg2 += t_i[np.where(u_i == ticker_dict[i])][-1]/total_events
print('')
print(f'Total number of total events for 10 companies: ',avg1)

print(f'Average number of total events for 10 companies: ',avg1/10)
print(f'Average time between each event for 10 companies: ',avg2/10)<jupyter_output>number of events of stock AAPL: 458
average time between each event of 'type' AAPL: 9.558951965065502
number of events of stock MSFT: 455
average time between each event of 'type' MSFT: 9.6
number of events of stock JPM: 453
average time between each event of 'type' JPM: 9.664459161147903
number of events of stock GS: 454
average time between each event of 'type' GS: 9.62114537444934
number of events of stock PFE: 453
average time between each event of 'type' PFE: 9.64459161147903
number of events of stock MRK: 453
average time between each event of 'type' MRK: 9.631346578366445
number of events of stock NKE: 455
average time between each event of 'type' NKE: 9.621978021978022
number of events of stock HD: 453
average time between each event of 'type' HD: 9.631346578366445
number of events of stock GOOGL: 453
average time between each event of 'type' GOOGL: 9.660044150110375
number of events of stock AMZN: 453
average time between each event of 'type' AMZN: 9.664459161147903

Total num[...]<jupyter_text>Collect relevant data such as event times (largest 20% price drop, only include days were return was below 20% quantile), markers,...<jupyter_code># Collect event times

t_i=[]
u_i=[] # collect marker for each event times

ticker_id = [] # collect integer marker of each stock, index correspodns to stock AABA, 1 correponds to stock AAPL
ticker = []
ticker_dict = {} # collect key-value pairs where key represents tikcer of a stock and it's corresponding value an integer marker


for u,col in enumerate(df[ten_companies]):
    Tdiff = df[col].diff()
    timestamps=dates[Tdiff<=Tdiff.quantile(0.2)] # return timestamps where subsequent price difference was less than 10% quantile
    t_i.extend((timestamps - pd.Timestamp(2006,1,3)).dt.days.astype(float)) # measure time 
    u_i.extend(np.repeat(u,len(timestamps)))
    ticker_dict[col] = u
    ticker_id += [u]
    ticker += [col]

t_i = np.array(t_i)
#t_i=np.array(t_i)+np.random.rand(len(t_i))
u_i=np.array(u_i)
perm = np.argsort(t_i)
t_i = t_i[perm] # superposed timestamps
u_i = u_i[perm]
tstamps_ten_univ = []
tstamps_ten_univ_dates = []

count = 0
for i in ticker:
    tdiff_i = df[i].diff()
    tstamps_ten_univ_dates += [dates[tdiff_i<=tdiff_i.quantile(0.2)]]
    tstamps_i = (tstamps_ten_univ_dates[count] - pd.Timestamp(2006,1,3)).dt.days.astype(float) # change time units from days to months
    tstamps_ten_univ += [tstamps_i]
    count += 1
ts_i_with_dates = []
for i in range(len(ten_companies)):
    ts_i_with_dates += [pd.DataFrame({'Date': tstamps_ten_univ_dates[i].values, 'event time': tstamps_ten_univ[i].values})]
    #ts_i_with_dates[i].to_csv('ten_stocks_event_times/'+ten_companies[i]+'_ts_with_dates.csv', index=None)
print(ts_i_with_dates[0]['event time'])
avg1 = 0
avg2 = 0

for i in ticker:
    total_events = t_i[np.where(u_i == ticker_dict[i])].shape[0]
    print(f'number of events of stock '+i+':', total_events)
    print(f'average time between each event of \'type\' '+i+':', t_i[np.where(u_i == ticker_dict[i])][-1]/total_events)
    avg1 += total_events
    avg2 += t_i[np.where(u_i == ticker_dict[i])][-1]/total_events
print('')

print(f'Total number of total events for 10 companies: ',avg1)

print(f'Average number of total events for 10 companies: ',avg1/10)
print(f'Average time between each event for 10 companies: ',avg2/10)<jupyter_output>number of events of stock AAPL: 613
average time between each event of 'type' AAPL: 7.141924959216966
number of events of stock MSFT: 604
average time between each event of 'type' MSFT: 7.231788079470198
number of events of stock JPM: 613
average time between each event of 'type' JPM: 7.141924959216966
number of events of stock GS: 605
average time between each event of 'type' GS: 7.224793388429752
number of events of stock PFE: 610
average time between each event of 'type' PFE: 7.163934426229508
number of events of stock MRK: 604
average time between each event of 'type' MRK: 7.223509933774834
number of events of stock NKE: 611
average time between each event of 'type' NKE: 7.1653027823240585
number of events of stock HD: 606
average time between each event of 'type' HD: 7.1996699669967
number of events of stock GOOGL: 605
average time between each event of 'type' GOOGL: 7.234710743801653
number of events of stock AMZN: 604
average time between each event of 'type' AMZN: 7.24834437086[...]<jupyter_text>I havent saved data where time is measured in months instead of days<jupyter_code>pd.DataFrame(data, columns=['event time','node']).to_csv('multivariate_ts_data.csv', index=None)
node = 0
ts_col = [np.array([])]
for i in range(len(ticker)):
    ts_col = np.append(ts_col, data[np.where(data[:,1] == node)][:,0])
    node += 1

node = 0
ts_col = []
for i in range(len(ticker)):
    ts_col += [data[np.where(data[:,1] == node)][:,0].flatten()]
    node += 1
ts_col = np.array(ts_col)
print(ticker)
pd.DataFrame(data[np.where(data[:,1] == 4)][:,0], columns=['event time']).to_csv('Pfizer_ts_data.csv', index=False)
HP.plot_event_times(data, num_of_nodes=10,  Time_horizon=-999, company_ticker=ticker, xaxislabel=r't, number of days since January 3, 2006', show_time_periods=False, labeled=True)
# plt.show()
plt.ion()<jupyter_output><empty_output><jupyter_text>20% drop<jupyter_code># Collect event times

t_i20=[]
u_i20=[] # collect marker for each event times

ticker_id = [] # collect integer marker of each stock, index correspodns to stock AABA, 1 correponds to stock AAPL
ticker = []
ticker_dict = {} # collect key-value pairs where key represents tikcer of a stock and it's corresponding value an integer marker


for u,col in enumerate(df[ten_companies]):
    Tdiff = np.log(df[col]).diff()
    timestamps=dates[Tdiff<=Tdiff.quantile(0.2)] # return timestamps where subsequent price difference was less than 0.2
    t_i20.extend((timestamps - pd.Timestamp(2006,1,3)).dt.days.astype(float))
    u_i20.extend(np.repeat(u,len(timestamps)))
    ticker_dict[col] = u
    ticker_id += [u]
    ticker += [col]


t_i20=np.array(t_i20)+np.random.rand(len(t_i20))
u_i20=np.array(u_i20)
perm = np.argsort(t_i20)
t_i20 = t_i20[perm] # superposed timestamps
u_i20 = u_i20[perm]



for i in ticker:
    total_events = t_i20[np.where(u_i20 == ticker_dict[i])].shape[0]
    print(f'number of events of stock '+i+':', total_events)
    avg = t_i20[np.where(u_i20 == ticker_dict[i])][-1]/total_events
    print(f'average time between each event of \'type\' '+i+':')
    print(r'%.3f days' % avg)<jupyter_output>number of events of stock AAPL: 604
average time between each event of 'type' AAPL:
7.245 days
number of events of stock MSFT: 604
average time between each event of 'type' MSFT:
7.207 days
number of events of stock JPM: 604
average time between each event of 'type' JPM:
7.223 days
number of events of stock GS: 604
average time between each event of 'type' GS:
7.233 days
number of events of stock PFE: 604
average time between each event of 'type' PFE:
7.234 days
number of events of stock MRK: 604
average time between each event of 'type' MRK:
7.225 days
number of events of stock NKE: 604
average time between each event of 'type' NKE:
7.246 days
number of events of stock HD: 604
average time between each event of 'type' HD:
7.211 days
number of events of stock GOOGL: 604
average time between each event of 'type' GOOGL:
7.208 days
number of events of stock AMZN: 604
average time between each event of 'type' AMZN:
7.249 days<jupyter_text>90% price jump<jupyter_code># Collect event times

t_i90=[]
u_i90=[] # collect marker for each event times

ticker_id = [] # collect integer marker of each stock, index correspodns to stock AABA, 1 correponds to stock AAPL
ticker = []
ticker_dict = {} # collect key-value pairs where key represents tikcer of a stock and it's corresponding value an integer marker


for u,col in enumerate(df[ten_companies]):
    Tdiff = (df[col]).diff()
    timestamps=dates[Tdiff>=Tdiff.quantile(0.9)] # return timestamps where subsequent price difference was less than 0.2
    t_i90.extend((timestamps - pd.Timestamp(2006,1,3)).dt.days.astype(float))
    u_i90.extend(np.repeat(u,len(timestamps)))
    ticker_dict[col] = u
    ticker_id += [u]
    ticker += [col]


t_i90=np.array(t_i90)+np.random.rand(len(t_i90))
u_i90=np.array(u_i90)
perm = np.argsort(t_i90)
t_i90 = t_i90[perm] # superposed timestamps
u_i90 = u_i90[perm]

for i in ticker:
    total_events = t_i90[np.where(u_i90 == ticker_dict[i])].shape[0]
    print(f'number of events of stock '+i+':', total_events)
    print(f'average time between each event of \'type\' '+i+':', t_i90[np.where(u_i90 == ticker_dict[i])][-1]/total_events)<jupyter_output>number of events of stock AAPL: 303
average time between each event of 'type' AAPL: 14.414080194724823
number of events of stock MSFT: 303
average time between each event of 'type' MSFT: 14.403819918473534
number of events of stock JPM: 304
average time between each event of 'type' JPM: 14.375039087571116
number of events of stock GS: 302
average time between each event of 'type' GS: 14.471115837880753
number of events of stock PFE: 306
average time between each event of 'type' PFE: 14.262688718931082
number of events of stock MRK: 303
average time between each event of 'type' MRK: 14.39595108595655
number of events of stock NKE: 304
average time between each event of 'type' NKE: 14.376841098160119
number of events of stock HD: 302
average time between each event of 'type' HD: 14.489696643966676
number of events of stock GOOGL: 302
average time between each event of 'type' GOOGL: 14.461547275605175
number of events of stock AMZN: 302
average time between each event of 'type' AMZN: 14.4[...]<jupyter_text>Use Kolmogorov-Smirnov test to check if the 10 univariates follow homogenous Poisson process<jupyter_code>for u,col in enumerate(df[ten_companies]):
    Tdiff = df[col].diff()
    timestamps=dates[Tdiff<=Tdiff.quantile(0.1)] # return timestamps where subsequent price difference was less than 10% quantile
    t_i.extend((timestamps - pd.Timestamp(2006,1,3)).dt.days.astype(float)) # measure time 
    u_i.extend(np.repeat(u,len(timestamps)))
    ticker_dict[col] = u
    ticker_id += [u]
    ticker += [col]
w = np.linspace(1,12,100)
P = MHP.MHP() # instantiate MHP object
P.generate_seq(60) # mu=[0.1], alpha=[[0.5]], and omega=1.0

#P.data = data[np.where(data[:,1] == 0)]
#P.data[:,1] = 0
#print(data[np.where(data[:,1] == 3)])
P.mu = [0.1]
P.alpha = [0.3]
P.omega = 0.5


# P.plot_events()
# plt.show()

seed = 10

rng = np.random.default_rng(seed)
ahat = rng.uniform(0,1, size=(1,1))
mhat = rng.uniform(0,1, size=11)



#w = np.linspace(0.1,4.25,100)
ahat, mhat
ahat_arr = []; mhat_arr = []; LL_arr = []; res = []
u_i = np.linspace(0,0,np.shape(tstamps_PFE)[0])
data = []
counter = 0
for t,u in zip(tstamps_PFE, u_i):
    data.append([])
    data[counter] += [t, np.array([u], dtype='int')]
    counter += 1
data = np.array(data)
for w_hyperparam in range(len(w)):
    #print(f'hyperparameter omega: ', w[i])
    P.data = data[0:159,:]
    P.data[:,1] = 0
    res += [P.EM(np.array([ahat[0]]), np.array([mhat[0]]), w[w_hyperparam], verbose=False, seed1=None,seed2=None)]
    ahat_arr += [res[w_hyperparam][0]] 
    mhat_arr += [res[w_hyperparam][1]]
    LL_arr += [-res[w_hyperparam][-1]]
min_index = np.argmin(LL_arr)
omega_index = [min_index]
mhat_vals = mhat_arr[min_index]
ahat_vals = ahat_arr[min_index]
w[omega_index], mhat_vals, ahat_vals
plt.plot(w,LL_arr)
plt.show()
plt.ion()
seed = 100

rng = np.random.default_rng(seed)
ahat = rng.uniform(0,1, size=(1,1))
mhat = rng.uniform(0,1, size=1)



#w = np.linspace(0.1,4.25,100)
ahat, mhat
mu_val, alpha_tilde_val, beta_val, NLL = HP.EM(mhat, ahat, 3, data[:,0], Maxiter=10, num_of_tstamp = -1)
mu_val, alpha_tilde_val/beta_val, beta_val, NLL<jupyter_output><empty_output><jupyter_text>Pfizer data outlier investigation<jupyter_code>df2 = pd.read_csv('dataset-djia-2018-subset2.csv')
dates = pd.to_datetime(df2['Date']) # set dates as the Date of closing price column
df2[['PFE','MRK']].plot()
dates = pd.to_datetime(df2['Date']) # set dates as the Date of closing price column
df2.reset_index().plot(x='Date', y=['PFE','MRK'])
#plt.yscale('log')
t_i = []
u_i = []
df2 = df2.set_index('Date')
tdiff_PFE = df['PFE'].diff()
ts_PFE = dates[tdiff_PFE<=tdiff_PFE.quantile(0.1)]
tstamps_PFE = (ts_PFE - pd.Timestamp(2006,1,3)).dt.days.astype(float)
ts_PFE_date = pd.DataFrame({'Date': ts_PFE.values, 'event time': tstamps_PFE.values})<jupyter_output><empty_output><jupyter_text>Appled data investigation<jupyter_code>tdiff_AAPL = df['AAPL'].diff()
ts_AAPL = dates[tdiff_AAPL<=tdiff_AAPL.quantile(0.1)]
tstamps_AAPL = (ts_AAPL - pd.Timestamp(2006,1,3)).dt.days.astype(float)
AAPL_event_times = pd.DataFrame({'event time': tstamps_AAPL.values})
AAPL_event_times
AAPL_event_times.to_csv('AAPL_event_times.dat', header=None, index=None)<jupyter_output><empty_output><jupyter_text>Save event times with dates data of each stock<jupyter_code>tstamps_ten_univ = []
tstamps_ten_univ_dates = []

count = 0
for i in ticker:
    tdiff_i = df[i].diff()
    tstamps_ten_univ_dates += [dates[tdiff_i<=tdiff_i.quantile(0.1)]]
    tstamps_i = (tstamps_ten_univ_dates[count] - pd.Timestamp(2006,1,3)).dt.days.astype(float)
    tstamps_ten_univ += [tstamps_i]
    count += 1
ts_i_with_dates = []
for i in range(len(ticker)):
    ts_i_with_dates += [pd.DataFrame({'Date': tstamps_ten_univ_dates[i].values, 'event time': tstamps_ten_univ[i].values})]
    ts_i_with_dates[i].to_csv('ten_stocks_event_times/'+ticker[i]+'_ts_with_dates.csv', index=None)
ts_MRK.iloc[0:100]
tstamps_MRK.iloc[0:100]
ts_PFE.iloc[0:160]
tstamps_PFE.iloc[0:160]
tstamps_PFE_subset = np.array(tstamps_PFE, dtype=float)[0:159]
tstamps_PFE_subset
from scipy import stats
delta_t = stats.expon.rvs(size=10000)
t = np.cumsum(delta_t)
stats.kstest(t/t.max(), 'uniform')
plt.hist(t/t.max())
plt.show()
plt.ion()
delta_t = np.diff(tstamps_PFE_check) # waiting times between events
#delta_t = np.diff(data[:,0])
t = np.cumsum(delta_t) # cumulative sum of the waiting times
sample = t/t.max()
sample = np.array(sample, dtype=float)
stats.kstest(sample, 'uniform')
plt.hist(sample)
stats_res = []
counter = 0
print("Kolmogorov-Smirnov test to check if the 10 univariates follow homogenous Poisson process")

for i in ticker:
    delta_t = np.diff(t_i[np.where(u_i == ticker_dict[i])])
    t = np.cumsum(delta_t) # cumulative sum of the waiting times
    sample = t/t.max()
    sample = np.array(sample, dtype=float)
    stats_res += [stats.kstest(sample, 'uniform')]
    print("")
    print("K-S test result of "+i+':')
    print(f'p-value =',stats_res[counter][1])
    counter += 1
for i in ticker:
    total_events = t_i20[np.where(u_i20 == ticker_dict[i])].shape[0]
    print(f'number of events of stock '+i+':', total_events)
    avg = t_i20[np.where(u_i20 == ticker_dict[i])][-1]/total_events
    print(f'average time between each event of \'type\' '+i+':')
    print(r'%.3f days' % avg)
t_i[np.where(u_i == ticker_dict[i])].shape<jupyter_output><empty_output><jupyter_text>Use EM to estimate alpha and mu given timestamps of univariate HP that models the events of each company using initial guess of alpha and mu   Compute log-likelihoods of each selected stock by tuning hyperparameter omega to find omega that gives the highest value of log-likelihood using EM<jupyter_code>P = MHP.MHP() # instantiate MHP object
P.generate_seq(60) # mu=[0.1], alpha=[[0.5]], and omega=1.0

#P.data = data[np.where(data[:,1] == 0)]
#P.data[:,1] = 0
#print(data[np.where(data[:,1] == 3)])
P.mu = [0.1]
P.alpha = [0.3]
P.omega = 0.5


# P.plot_events()
# plt.show()

seed = 1

rng = np.random.default_rng(seed)
ahat = rng.uniform(0,1, size=(10,10))
mhat = rng.uniform(0,1, size=10)



#w = np.linspace(0.1,4.25,100)
ahat, mhat
w = np.linspace(0.05,4,100) # hyperparameter values chosen for sufficient coverage
w<jupyter_output><empty_output><jupyter_text>Hyperparameter optimisation  Vary hyperparmaeter $\omega$ to compute log-likelihoods using EM algorithm and find the "best" parametric estimate of $\omega$, i.e., $\omega$ with least value of log-likelihood<jupyter_code>omega_index = [] # store index of best value of omega 
ahat_vals = [] 
mhat_vals = []

ahat_arr = []; mhat_arr = []
LL_arr = []
# axes are in a two-dimensional array, indexed by [row, col]
fig, axs = plt.subplots(2, 5)

node = 0
for i in range(2):
    axs[i, 0].set_ylabel('Negative log-likelihood')

    for j in range(5):
        
        ahat_arr += [[]]; mhat_arr += [[]]
        LL_arr += [[]]; res = []
        for w_hyperparam in range(len(w)):
            #print(f'hyperparameter omega: ', w[i])
            P.data = data[np.where(data[:,1] == node)]
            P.data[:,1] = 0
            res += [P.EM(np.array([ahat[node]]), np.array([mhat[node]]), w[w_hyperparam], verbose=False, seed1=99,seed2=100)]
            ahat_arr[node] += [res[w_hyperparam][0]] 
            mhat_arr[node] += [res[w_hyperparam][1]]
            LL_arr[node] += [-res[w_hyperparam][-1]]
            #print(estimates[i])
            #print('\n')
        
        
        
#         while abs(LL_arr[count] - LL_arr[-1]) >= epsilon:
#             count += 1
#         print(abs(LL_arr[count] - LL_arr[-1]))
        min_index = np.argmin(LL_arr[node])
        omega_index += [min_index]
        mhat_vals += [mhat_arr[node][min_index]]
        ahat_vals += [ahat_arr[node][min_index]]
        
        ax = axs[i, j]
        ax.plot(w, LL_arr[node])
        ax.text(0.7, 0.25, ten_companies[node],transform=ax.transAxes,
            fontsize=18,verticalalignment ='bottom',
            horizontalalignment ='right', color ='green',)
        #axs[i, j].axvline(x=w[count], ymin=LL_arr[count], ymax=LL_arr[0], label=r'$\hat{\omega} = $'+str(w[count]), c='b', linestyle='dashed')
        
        ax.scatter(x=w[min_index],y=LL_arr[node][min_index], c='b',s=80,marker='X',linewidths=1.5)
        #ax.plot(w[count],LL_arr[count],color='blue',label=r'$\hat{\omega} = $'+str(w[count]),linewidth=2, markersize=50)

        node += 1
        print('Working')

        #ax[i, j].text(0.5, 0.5, str((i, j)), fontsize=18, ha='center')
for j in range(5):
    axs[1, j].set_xlabel(r'omega, $\omega$')

plt.tight_layout()
#plt.savefig('LL_fig_iter6.png')
plt.show()
plt.ion()
omega_index
for i in mhat_vals:
    print(r'mu: %.3f' % i[0])
for i in ahat_vals:
    print(r'alpha: %.3f' % i[0])
with np.printoptions(precision=3, suppress=True):
    print(w[omega_index])
squad
sqad<jupyter_output>[array([0.0692]), array([0.0698]), array([0.0692]), array([0.0694]), array([0.0693]), array([0.0697]), array([0.069]), array([0.0693]), array([0.0694]), array([0.069])]<jupyter_text>Estimated parameter values:<jupyter_code>omega_index, omega_index1
with np.printoptions(precision=3, suppress=True):
    print(omega_index)

ahat_vals, ahat_vals1
with np.printoptions(precision=4, suppress=True):
    print(mhat_vals)
    print(mhat_vals1)
ahat_vals, ahat_vals1
B
HP.plot_event_times(data, num_of_nodes=10,  Time_horizon=-999, company_ticker=ticker, xaxislabel=r't, number of days since January 3, 2006', show_time_periods=False, labeled=True)
seed = 99
ts = np.array([])
for i in range(2):
    np.append(ts, HP.simulate_timestamps_till_horizon(mhat_vals[i], ahat_vals[i], w[i] , Thorizon = np.amax(data[:,0]), seed=seed, node=i, output_rejected_data=False))
ts
    
ts1 = HP.simulate_timestamps_till_horizon(mhat_vals[0], ahat_vals[0], w[0] , Thorizon = np.amax(data[:,0]), seed=None, output_rejected_data=False)
ts2 = HP.simulate_timestamps_till_horizon(mhat_vals[1], ahat_vals[0], w[1] , Thorizon = np.amax(data[:,0]), seed=None, output_rejected_data=False)
ts3 = HP.simulate_timestamps_till_horizon(mhat_vals[2], ahat_vals[2], w[2] , Thorizon = np.amax(data[:,0]), seed=None, output_rejected_data=False)
ts4 = HP.simulate_timestamps_till_horizon(mhat_vals[3], ahat_vals[3], w[3] , Thorizon = np.amax(data[:,0]), seed=None, output_rejected_data=False)
ts5 = HP.simulate_timestamps_till_horizon(mhat_vals[4], ahat_vals[4], w[4] , Thorizon = np.amax(data[:,0]), seed=None, output_rejected_data=False)

#a_hat_1 = 0.00882517; mhat_vals_1 = 0.06859064; w_1 = 1.5848484848484852

plt.plot(ts1,[0]*len(ts1), marker='o')
plt.plot(ts2,[-1]*len(ts2), marker='o')
plt.plot(ts3,[-2]*len(ts3), marker='o')
plt.plot(ts4,[-3]*len(ts4), marker='o')
plt.plot(ts5,[-4]*len(ts5), marker='o')
B<jupyter_output><empty_output><jupyter_text>Save this session and load it back<jupyter_code>import dill
dill.dump_session('DJIA_univHP_env_September_version2.db')
#dill.load_session('DJIA_univHP_env_September_version2.db')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13958333333333334, 'fraction_numerical': 0.05438218390804598, 'mean_word_length': 4.952747487705794, 'pattern_counts': {'":': 0, '<': 74, '<?xml version=': 0, '>': 65, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '609436', 'license': 'MIT', 'n_tokens_mistral': 11672, 'n_tokens_neox': 10585, 'n_words': 2673, 'path': 'Financial_data_results/univariate_case/DJIA_univariate_latest.ipynb', 'repo_name': 'rai-d-babi/hawkesProcess', 'size': '27840'}