Simulation -- Birthday Paradox Market

import numpy as np ## Load Numpy import matplotlib.pyplot as plt ## Load matplotlib @interact def plot_birthday_paradox_market(K = slider(15,100,1,label='# Firms/Industry'), I = slider(15,50,1,label='# Industries'), delta = slider(0.001,0.05,0.001,default=0.05,label='Price Impact'), sig_z = slider(0.005,0.05,0.005,default=0.025,label='Noise Volatility'), N = slider(10,25,1,label='Search Depth') ): T = 1000 ## Number of periods I = int(np.round(I)) ## Number of industries K = int(np.round(K)) ## Number of assets per industry A = K * I ## Number of assets N = int(np.round(N)) ## Number of performing assets if (N > I): print "Warning!!! Guaranteed to find surprise: (N > I)!!!" v = np.zeros((T, A)) ## (T x A) matrix of excess asset returns s_up = np.zeros((T, I)) ## (T x I) matrix of number of positive surprises s_dn = np.zeros((T, I)) ## (T x I) matrix of number of negative surprises order_dict = {} ## Dict of order of previous period's returns ################################################################################################## ################################################################################################## ################################################################################################## v[0,:] = np.random.normal(0, sig_z, A) ## Initialize first period's returns for t in range(1, T): ## Loop over periods 1 to T ## Loop over all A assets at time t ## Fill up dictionary of returns ## Sort dictionary by key value for a in range(0,A): order_dict[a] = v[t-1,a] top_N = sorted(order_dict, key = order_dict.__getitem__, reverse = True)[0:N] bot_N = sorted(order_dict, key = order_dict.__getitem__, reverse = False)[0:N] ## Loop over performing assets ## Identify industry of asset ## Add 1 to surprise matrix for industry i at time t for n in range(0,N): for i in range(0,I): if (top_N[n] < K * (i+1)) and (top_N[n] >= K * i): s_up[t,i] = s_up[t,i] + 1 if (bot_N[n] < K * (i+1)) and (bot_N[n] >= K * i): s_dn[t,i] = s_dn[t,i] + 1 ## Loop over industries ## Fill in random returns ## If industry experiences surprise, augment mean return by delta for i in range(0,I): if (s_up[t,i] >= 2): v[t,(K*i):(K*(i+1))] = np.random.normal(delta, sig_z, K) elif (s_dn[t,i] >= 2): v[t,(K*i):(K*(i+1))] = np.random.normal(-delta, sig_z, K) else: v[t,(K*i):(K*(i+1))] = np.random.normal(0, sig_z, K) ################################################################################################## ################################################################################################## ################################################################################################## ## Define plotting environment plt.figure(figsize=(8, 24)) plt.subplots_adjust(hspace=0.2) sub_title_1 = "Avg. Excess Returns by Industry" plt.subplot(611, title = sub_title_1) ## Create time series plot of excess returns for i in range(0, I): ## Loop over industries r_imean = [] for t in range(0, T): ## Loop over time periods r_imean.append(np.mean(v[t,(K*i):(K*(i+1))])) ## Compute industry i's mean return at time t plt.plot(range(0, T), r_imean) ## Plot industry i's mean from from time 0 to T plt.grid(True) sub_title_2 = "Avg. Excess Returns by Industry" plt.subplot(612, title = sub_title_2) ## Create histogram subplot r_imean = [] for i in range(0,I): ## Loop over industries for t in range(0, T): ## Loop over times r_imean.append(np.mean(v[t,(K*i):(K*(i+1))])) ## Compute industry mean at time t plt.hist(r_imean, ## Plot histogram of industry means 150, facecolor = 'green', alpha = 0.50, align = 'mid' ) sub_title_3 = "# of Positive s-Way Surprises" plt.subplot(613, title = sub_title_3) ## Create time series plot of surprises for i in range(0, I): ## Loop over industries plt.plot(range(0, T), s_up[:,i]) ## Plot industry i's surprise count from from time 0 to T plt.grid(True) sub_title_4 = "# of Positive s-Way Surprises" plt.subplot(614, title = sub_title_4) ## Create histogram subplot s_long = [] for i in range(0,I): ## Loop over industries for t in range(0,T): ## Loop over periods if (s_up[t,i] != 0): s_long.append(s_up[t,i]) ## Add industry counts to list plt.hist(s_long, ## Plot histogram of s-way surprises 100, facecolor = 'green', alpha = 0.50, align = 'mid' ) sub_title_5 = "Positive Surprise History {s_t, s_{t+1},... | s_t >= 2}" plt.subplot(615, title = sub_title_5) ## Create plot of lifetimes for i in range(0,I): ## Loop over industries s_sts = [] during_surprise = False for t in range(0,T): ## Loop over periods if (s_up[t,i] >= 2): during_surprise = True else: if (s_sts != []): plt.plot(range(0, len(s_sts)), s_sts) s_sts = [] during_surprise = False if (during_surprise == True): s_sts.append(s_up[t,i]) plt.grid(True) sub_title_6 = "Cumulative Hazard Rate Positive s-Way Surprise" plt.subplot(616, title = sub_title_6) ## Create plot of lifetimes survival_time = [] for i in range(0,I): ## Loop over industries s_sts = [] during_surprise = False for t in range(0,T): ## Loop over periods if (s_up[t,i] >= 2): during_surprise = True else: if (s_sts != []): survival_time.append(len(s_sts)) s_sts = [] during_surprise = False if (during_surprise == True): s_sts.append(s_up[t,i]) plt.hist(survival_time, ## Plot histogram of s-way surprises 100, facecolor = 'green', alpha = 0.50, cumulative = True, normed = True, align = 'mid' ) plt.savefig('plt.png') plt.clf() ## Update plot.