diff --git a/.gitignore b/.gitignore index 3eb56cb..d21a2d5 100644 --- a/.gitignore +++ b/.gitignore @@ -150,3 +150,9 @@ cython_debug/ # and can be added to the global gitignore or merged into this file. For a more nuclear # option (not recommended) you can uncomment the following to ignore the entire idea folder. .idea/ + +# PNG images +*.png + +# Mac stuff +*.DS_Store \ No newline at end of file diff --git a/7.26 comp 1000/initial window and V_hat_comp_small0.01.png b/7.26 comp 1000/initial window and V_hat_comp_small0.01.png deleted file mode 100644 index fc043a1..0000000 Binary files a/7.26 comp 1000/initial window and V_hat_comp_small0.01.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and V_hat_comp_small0.02.png b/7.26 comp 1000/initial window and V_hat_comp_small0.02.png deleted file mode 100644 index fc043a1..0000000 Binary files a/7.26 comp 1000/initial window and V_hat_comp_small0.02.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and V_hat_comp_small0.03.png b/7.26 comp 1000/initial window and V_hat_comp_small0.03.png deleted file mode 100644 index fc043a1..0000000 Binary files a/7.26 comp 1000/initial window and V_hat_comp_small0.03.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and age of participants_comp_small0.01.png b/7.26 comp 1000/initial window and age of participants_comp_small0.01.png deleted file mode 100644 index 475e724..0000000 Binary files a/7.26 comp 1000/initial window and age of participants_comp_small0.01.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and age of participants_comp_small0.02.png b/7.26 comp 1000/initial window and age of participants_comp_small0.02.png deleted file mode 100644 index 07fdb9d..0000000 Binary files a/7.26 comp 1000/initial window and age of participants_comp_small0.02.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and age of participants_comp_small0.03.png b/7.26 comp 1000/initial window and age of participants_comp_small0.03.png deleted file mode 100644 index 671a605..0000000 Binary files a/7.26 comp 1000/initial window and age of participants_comp_small0.03.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and cst belief of participants in age groups_comp_small0.01.png b/7.26 comp 1000/initial window and cst belief of participants in age groups_comp_small0.01.png deleted file mode 100644 index c84177a..0000000 Binary files a/7.26 comp 1000/initial window and cst belief of participants in age groups_comp_small0.01.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and cst belief of participants in age groups_comp_small0.02.png b/7.26 comp 1000/initial window and cst belief of participants in age groups_comp_small0.02.png deleted file mode 100644 index 06c4883..0000000 Binary files a/7.26 comp 1000/initial window and cst belief of participants in age groups_comp_small0.02.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and cst belief of participants in age groups_comp_small0.03.png b/7.26 comp 1000/initial window and cst belief of participants in age groups_comp_small0.03.png deleted file mode 100644 index 429eb4e..0000000 Binary files a/7.26 comp 1000/initial window and cst belief of participants in age groups_comp_small0.03.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and cutoff belief_comp_small0.01.png b/7.26 comp 1000/initial window and cutoff belief_comp_small0.01.png deleted file mode 100644 index ee8fd60..0000000 Binary files a/7.26 comp 1000/initial window and cutoff belief_comp_small0.01.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and cutoff belief_comp_small0.02.png b/7.26 comp 1000/initial window and cutoff belief_comp_small0.02.png deleted file mode 100644 index 5e242dd..0000000 Binary files a/7.26 comp 1000/initial window and cutoff belief_comp_small0.02.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and cutoff belief_comp_small0.03.png b/7.26 comp 1000/initial window and cutoff belief_comp_small0.03.png deleted file mode 100644 index df9e313..0000000 Binary files a/7.26 comp 1000/initial window and cutoff belief_comp_small0.03.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and interest rate_comp_small0.01.png b/7.26 comp 1000/initial window and interest rate_comp_small0.01.png deleted file mode 100644 index 2136708..0000000 Binary files a/7.26 comp 1000/initial window and interest rate_comp_small0.01.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and interest rate_comp_small0.02.png b/7.26 comp 1000/initial window and interest rate_comp_small0.02.png deleted file mode 100644 index 607e19c..0000000 Binary files a/7.26 comp 1000/initial window and interest rate_comp_small0.02.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and interest rate_comp_small0.03.png b/7.26 comp 1000/initial window and interest rate_comp_small0.03.png deleted file mode 100644 index 9d694a7..0000000 Binary files a/7.26 comp 1000/initial window and interest rate_comp_small0.03.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and market price of risk_comp_small0.01.png b/7.26 comp 1000/initial window and market price of risk_comp_small0.01.png deleted file mode 100644 index 1d517b6..0000000 Binary files a/7.26 comp 1000/initial window and market price of risk_comp_small0.01.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and market price of risk_comp_small0.02.png b/7.26 comp 1000/initial window and market price of risk_comp_small0.02.png deleted file mode 100644 index ebddf0f..0000000 Binary files a/7.26 comp 1000/initial window and market price of risk_comp_small0.02.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and market price of risk_comp_small0.03.png b/7.26 comp 1000/initial window and market price of risk_comp_small0.03.png deleted file mode 100644 index ce9044d..0000000 Binary files a/7.26 comp 1000/initial window and market price of risk_comp_small0.03.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and number of cohorts_comp_small0.01.png b/7.26 comp 1000/initial window and number of cohorts_comp_small0.01.png deleted file mode 100644 index 4e2ef4d..0000000 Binary files a/7.26 comp 1000/initial window and number of cohorts_comp_small0.01.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and number of cohorts_comp_small0.02.png b/7.26 comp 1000/initial window and number of cohorts_comp_small0.02.png deleted file mode 100644 index e8c3608..0000000 Binary files a/7.26 comp 1000/initial window and number of cohorts_comp_small0.02.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and number of cohorts_comp_small0.03.png b/7.26 comp 1000/initial window and number of cohorts_comp_small0.03.png deleted file mode 100644 index f54c014..0000000 Binary files a/7.26 comp 1000/initial window and number of cohorts_comp_small0.03.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and participation rate_comp_small0.01.png b/7.26 comp 1000/initial window and participation rate_comp_small0.01.png deleted file mode 100644 index 03961cb..0000000 Binary files a/7.26 comp 1000/initial window and participation rate_comp_small0.01.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and participation rate_comp_small0.02.png b/7.26 comp 1000/initial window and participation rate_comp_small0.02.png deleted file mode 100644 index 11b734c..0000000 Binary files a/7.26 comp 1000/initial window and participation rate_comp_small0.02.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and participation rate_comp_small0.03.png b/7.26 comp 1000/initial window and participation rate_comp_small0.03.png deleted file mode 100644 index 208433d..0000000 Binary files a/7.26 comp 1000/initial window and participation rate_comp_small0.03.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and wealth share in age groups_comp_small0.01.png b/7.26 comp 1000/initial window and wealth share in age groups_comp_small0.01.png deleted file mode 100644 index e057cda..0000000 Binary files a/7.26 comp 1000/initial window and wealth share in age groups_comp_small0.01.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and wealth share in age groups_comp_small0.02.png b/7.26 comp 1000/initial window and wealth share in age groups_comp_small0.02.png deleted file mode 100644 index 3590722..0000000 Binary files a/7.26 comp 1000/initial window and wealth share in age groups_comp_small0.02.png and /dev/null differ diff --git a/7.26 comp 1000/initial window and wealth share in age groups_comp_small0.03.png b/7.26 comp 1000/initial window and wealth share in age groups_comp_small0.03.png deleted file mode 100644 index 7925ddb..0000000 Binary files a/7.26 comp 1000/initial window and wealth share in age groups_comp_small0.03.png and /dev/null differ diff --git a/V_hat_experiment.py b/V_hat_experiment.py index cbcfb48..05f58c4 100644 --- a/V_hat_experiment.py +++ b/V_hat_experiment.py @@ -7,13 +7,13 @@ # todo: the connection between belief and wealth? # Learning from repeated negative economic shocks: lead to both worse wealth condition and pessimism # modes = ['drop', 'ric_free'] -mode = 'rich_free' +mode = "rich_free" # mode = 'drop' # mode = 'comp' # mode = 'keep' # zoom_in = 'small' # zoom_in = 'large' -zoom_in = 'overall' +zoom_in = "overall" a = np.arange(1, 13, 1) b = np.arange(24, 241, 12) Npres = np.append(a, b) @@ -66,8 +66,8 @@ dR_matrix = np.zeros((T_hat_dimension, nu_dimension, Mpaths)) -dZ_matrix = np.load('dZ_matrix.npy') -dZ_build_matrix = np.load('dZ_build_matrix.npy') +dZ_matrix = np.load("dZ_matrix.npy") +dZ_build_matrix = np.load("dZ_build_matrix.npy") # The main loop builds up the economy with a large number of cohorts, and simulates the stationary economy forward for l in range(Mpaths): @@ -76,7 +76,7 @@ for k, T_hat in enumerate(T_hats): Npre = int(Npres[k]) - Vhat = (sigma_Y ** 2) / T_hat # prior variance + Vhat = (sigma_Y**2) / T_hat # prior variance print(T_hat, Npre, Vhat) for m, nu in enumerate(nus): @@ -84,7 +84,9 @@ beta = rho + nu - tax tau = np.arange(T_cohort, 0, -dt) # age from 500 to 0 - cohort_size = nu * np.exp(-nu * (tau - dt)) * dt # cohort size when a new cohort is just born + cohort_size = ( + nu * np.exp(-nu * (tau - dt)) * dt + ) # cohort size when a new cohort is just born # create age quartiles for analysis cummu_popu = np.cumsum(cohort_size) @@ -93,7 +95,7 @@ tau_cutoff3 = np.searchsorted(cummu_popu, 0.25) cutoffs = [Nc, tau_cutoff1, tau_cutoff2, tau_cutoff3, 0] - if mode == 'drop' or mode == 'keep' or mode == 'comp': + if mode == "drop" or mode == "keep" or mode == "comp": ( r, theta, @@ -108,12 +110,30 @@ w_cohort, age_parti, n_parti, - ) = simulate(mode, Nc, Nt, dt, rho, nu, Vhat, mu_Y, sigma_Y, sigma_S, tax, beta, Npre, Ninit, T_hat, - dZ_build, dZ, tau, - cohort_size) + ) = simulate( + mode, + Nc, + Nt, + dt, + rho, + nu, + Vhat, + mu_Y, + sigma_Y, + sigma_S, + tax, + beta, + Npre, + Ninit, + T_hat, + dZ_build, + dZ, + tau, + cohort_size, + ) invest_tracker = pi > 0 - elif mode == 'rich_free' or mode == 'back_collect' or mode == 'back_renew': + elif mode == "rich_free" or mode == "back_collect" or mode == "back_renew": ( r, theta, @@ -141,10 +161,29 @@ Delta_bar_parti, Delta_bar_long, Delta_bar_short, - ) = simulate_partial_constraint(mode, Nc, Nt, dt, rho, nu, Vhat, mu_Y, sigma_Y, sigma_S, tax, beta, Npre, - Ninit, T_hat, dZ_build, dZ, tau, cohort_size) + ) = simulate_partial_constraint( + mode, + Nc, + Nt, + dt, + rho, + nu, + Vhat, + mu_Y, + sigma_Y, + sigma_S, + tax, + beta, + Npre, + Ninit, + T_hat, + dZ_build, + dZ, + tau, + cohort_size, + ) else: - print('Error! Mode not defined') + print("Error! Mode not defined") break dR_matrix[k, m, l] = np.mean(dR) r_matrix[k, m, l] = np.mean(r) @@ -154,11 +193,13 @@ f_parti_matrix[k, m, l] = np.mean(f_parti) parti_rate = invest_tracker * cohort_size - belief = (Delta * sigma_Y + mu_Y) + belief = Delta * sigma_Y + mu_Y belief_weights = f * dt for i in range(4): - popu_age_matrix[k, m, l, i] = np.mean(np.sum(parti_rate[:, cutoffs[i + 1]:], axis=1)) + popu_age_matrix[k, m, l, i] = np.mean( + np.sum(parti_rate[:, cutoffs[i + 1] :], axis=1) + ) # weights_zero = (np.sum(invest_tracker[:, cutoffs[i + 1]:cutoffs[i]], # axis=1) == 0) # no one from the age group is participating @@ -174,11 +215,11 @@ # weights=belief_weights[:, cutoffs[i + 1]:cutoffs[i]], axis=1) # ) belief_age_matrix[k, m, l, i] = np.mean( - belief[:, cutoffs[i + 1]:cutoffs[i]] + belief[:, cutoffs[i + 1] : cutoffs[i]] ) wealthshare_age_matrix[k, m, l, i] = np.mean( - np.sum(f[:, cutoffs[i + 1]:cutoffs[i]] * dt, axis=1) + np.sum(f[:, cutoffs[i + 1] : cutoffs[i]] * dt, axis=1) ) age_parti_matrix[k, m, l] = np.mean(age_parti) @@ -189,26 +230,37 @@ # graphs: x = Npres -y0 = (np.ones(len(Npres)) * sigma_Y ** 2) / x +y0 = (np.ones(len(Npres)) * sigma_Y**2) / x y1 = np.mean(r_matrix, axis=2) y2 = np.mean(theta_matrix, axis=2) y3 = np.mean(popu_age_matrix, axis=2) y4 = np.mean(age_parti_matrix, axis=2) y5 = np.mean(n_parti_matrix, axis=2) y6 = -y2 * sigma_Y + mu_Y -y7 = np.nanmean(belief_age_matrix, axis=2) # consumption-weighted beliefs for participants from each age group +y7 = np.nanmean( + belief_age_matrix, axis=2 +) # consumption-weighted beliefs for participants from each age group y8 = np.nanmean(wealthshare_age_matrix, axis=2) # wealth share each age group y9 = np.mean(f_parti_matrix, axis=2) y10 = np.mean(Delta_bar_parti_matrix, axis=2) -xlabels = ['V_hat', 'interest rate', 'market price of risk', 'participation rate', 'average age of participants', - 'number of cohorts', 'cutoff belief to participate', 'average belief in age groups', - 'wealth share in age groups', 'consumption share of participants', 'estimation error of participants'] +xlabels = [ + "V_hat", + "interest rate", + "market price of risk", + "participation rate", + "average age of participants", + "number of cohorts", + "cutoff belief to participate", + "average belief in age groups", + "wealth share in age groups", + "consumption share of participants", + "estimation error of participants", +] ys = [y0, y1, y2, y3, y4, y5, y6, y7, y8, y9, y10] - for i in range(len(ys)): for j in range(nu_dimension): nu = nus[j] @@ -219,27 +271,78 @@ y_nu = ys[i] y = y_nu[:, j] if i == 3: - ax.fill_between(x, y[:, 0], color='steelblue', linewidth=0.4, label='20 < Age <= 35, youngest quartile') - ax.fill_between(x, y[:, 1], y[:, 0], color='darkseagreen', linewidth=0.4, label='35 < Age <= 55') - ax.fill_between(x, y[:, 2], y[:, 1], color='moccasin', linewidth=0.4, label='55 < Age <= 89') - ax.fill_between(x, y[:, 3], y[:, 2], color='pink', linewidth=0.4, label='Age > 89, oldest quartile') + ax.fill_between( + x, + y[:, 0], + color="steelblue", + linewidth=0.4, + label="20 < Age <= 35, youngest quartile", + ) + ax.fill_between( + x, + y[:, 1], + y[:, 0], + color="darkseagreen", + linewidth=0.4, + label="35 < Age <= 55", + ) + ax.fill_between( + x, + y[:, 2], + y[:, 1], + color="moccasin", + linewidth=0.4, + label="55 < Age <= 89", + ) + ax.fill_between( + x, + y[:, 3], + y[:, 2], + color="pink", + linewidth=0.4, + label="Age > 89, oldest quartile", + ) plt.legend() elif i == 7 or i == 8: - ax.plot(x, y[:, 0], color='steelblue', linewidth=0.4, label='20 < Age <= 35, youngest quartile') - ax.plot(x, y[:, 1], color='darkseagreen', linewidth=0.4, label='35 < Age <= 55') - ax.plot(x, y[:, 2], color='moccasin', linewidth=0.4, label='55 < Age <= 89') - ax.plot(x, y[:, 3], color='pink', linewidth=0.4, label='Age > 89, oldest quartile') + ax.plot( + x, + y[:, 0], + color="steelblue", + linewidth=0.4, + label="20 < Age <= 35, youngest quartile", + ) + ax.plot( + x, y[:, 1], color="darkseagreen", linewidth=0.4, label="35 < Age <= 55" + ) + ax.plot(x, y[:, 2], color="moccasin", linewidth=0.4, label="55 < Age <= 89") + ax.plot( + x, + y[:, 3], + color="pink", + linewidth=0.4, + label="Age > 89, oldest quartile", + ) plt.legend() else: ax.plot(x, y) - ax.set_xlabel('initial window (months)') + ax.set_xlabel("initial window (months)") if i == 0 or i == 6: ax.set_ylabel(xlabels[i]) else: - ax.set_ylabel('mean ' + xlabels[i]) + ax.set_ylabel("mean " + xlabels[i]) - plt.savefig('initial window and ' + xlabels[i] + '_' + mode + '_' + zoom_in + str(nu) + '.png', dpi=200, format="png") + plt.savefig( + "initial window and " + + xlabels[i] + + "_" + + mode + + "_" + + zoom_in + + str(nu) + + ".png", + dpi=200, + format="png", + ) # plt.savefig('initial window and ' + xlabels[i] + '_' + mode + '.png', dpi=500, format="png") plt.show() plt.close() - diff --git a/dataframe.py b/dataframe.py index cd06439..1764e09 100644 --- a/dataframe.py +++ b/dataframe.py @@ -1,12 +1,12 @@ import os import pandas as pd -data_file_folder = '.\SOI Tax Stats' +data_file_folder = ".\SOI Tax Stats" df = [] for file in os.listdir(data_file_folder): - if file.endswith('.csv'): - print('Loading file {0}...'.format(file)) + if file.endswith(".csv"): + print("Loading file {0}...".format(file)) df.append(pd.read_csv(os.path.join(data_file_folder, file))) print(len(df)) @@ -14,4 +14,4 @@ df_master = pd.concat(df, axis=0) # df_master.to_csv('master file.csv', index = False) -n1mean = df_master['n1'].mean() \ No newline at end of file +n1mean = df_master["n1"].mean() diff --git a/graphs.py b/graphs.py index 37768ce..3b449f7 100644 --- a/graphs.py +++ b/graphs.py @@ -3,13 +3,13 @@ import numpy as np # define the colors: -color1 = 'black' -color2 = 'mediumblue' -color3 = 'darkgreen' -color4 = 'orange' -color5 = 'red' -color6 = 'b' -color7 = 'g' +color1 = "black" +color2 = "mediumblue" +color3 = "darkgreen" +color4 = "orange" +color5 = "red" +color6 = "b" +color7 = "g" ####################################### ########## ONE RANDOM PATH ############ @@ -39,19 +39,19 @@ y13 = Delta_time_series[2] fig, ax1 = plt.subplots(figsize=(10, 5)) -ax1.set_xlabel('Time in simulation, one random path') -ax1.set_ylabel('Zt', color=color5) +ax1.set_xlabel("Time in simulation, one random path") +ax1.set_ylabel("Zt", color=color5) ax1.plot(t, y0, color=color5, linewidth=0.5) -ax1.tick_params(axis='y', labelcolor=color5) +ax1.tick_params(axis="y", labelcolor=color5) ax2 = ax1.twinx() -ax2.set_ylabel('Bias in belief and learning', color=color2) +ax2.set_ylabel("Bias in belief and learning", color=color2) ax2.set_ylim([-0.5, 0.5]) ax2.plot(t, y11, color=color2, linewidth=0.4) ax2.plot(t, y12, color=color3, linewidth=0.4) ax2.plot(t, y13, color=color4, linewidth=0.4) -ax2.tick_params(axis='y', labelcolor=color2) +ax2.tick_params(axis="y", labelcolor=color2) fig.tight_layout() # otherwise the right y-label is slightly clipped -plt.savefig('Zt and bias time series' + '.png', dpi=500) +plt.savefig("Zt and bias time series" + ".png", dpi=500) plt.show() ####################################### @@ -62,44 +62,60 @@ y22 = theta_drop fig, ax1 = plt.subplots(figsize=(15, 5)) -ax1.set_xlabel('Time in simulation, one random path') -ax1.set_ylabel('Zt', color=color5) +ax1.set_xlabel("Time in simulation, one random path") +ax1.set_ylabel("Zt", color=color5) ax1.plot(t, y0, color=color5, linewidth=0.5) -ax1.tick_params(axis='y', labelcolor=color5) +ax1.tick_params(axis="y", labelcolor=color5) ax2 = ax1.twinx() -ax2.set_ylabel('Market price of risk', color=color2) +ax2.set_ylabel("Market price of risk", color=color2) ax2.set_ylim([-1, 1]) -ax2.plot(t, y21, color=color2, linewidth=0.4, label='Complete market') -ax2.plot(t, y22, color=color3, linewidth=0.4, label='Short-sale constraint') -ax2.hlines(sigma_Y, xmin=0, xmax=500, color='purple', linestyles='--', linewidth=0.8, label='Representative agent') -ax2.tick_params(axis='y', labelcolor=color2) +ax2.plot(t, y21, color=color2, linewidth=0.4, label="Complete market") +ax2.plot(t, y22, color=color3, linewidth=0.4, label="Short-sale constraint") +ax2.hlines( + sigma_Y, + xmin=0, + xmax=500, + color="purple", + linestyles="--", + linewidth=0.8, + label="Representative agent", +) +ax2.tick_params(axis="y", labelcolor=color2) plt.legend() # fig.suptitle('Zt and Market Price of Risk') fig.tight_layout() # otherwise the right y-label is slightly clipped -plt.savefig('Zt and market price of risk' + '.png', dpi=500) +plt.savefig("Zt and market price of risk" + ".png", dpi=500) plt.show() # plot the interest rate y31 = r_comp y32 = r_drop -y33 = rho + mu_Y - sigma_Y ** 2 +y33 = rho + mu_Y - sigma_Y**2 fig, ax1 = plt.subplots(figsize=(15, 5)) -ax1.set_xlabel('Time in simulation, one random path') -ax1.set_ylabel('Zt', color=color5) +ax1.set_xlabel("Time in simulation, one random path") +ax1.set_ylabel("Zt", color=color5) ax1.plot(t, y0, color=color5, linewidth=0.5) -ax1.tick_params(axis='y', labelcolor=color5) +ax1.tick_params(axis="y", labelcolor=color5) ax2 = ax1.twinx() -ax2.set_ylabel('Interest rate (annual)', color=color2) +ax2.set_ylabel("Interest rate (annual)", color=color2) ax2.set_ylim([0, 0.05]) -ax2.plot(t, y31, color=color2, linewidth=0.4, label='Complete market') -ax2.plot(t, y32, color=color3, linewidth=0.4, label='Short-sale constraint') -ax2.hlines(y33, xmin=0, xmax=500, color='purple', linestyles='--', linewidth=0.8, label='Representative agent') -ax2.tick_params(axis='y', labelcolor=color2) +ax2.plot(t, y31, color=color2, linewidth=0.4, label="Complete market") +ax2.plot(t, y32, color=color3, linewidth=0.4, label="Short-sale constraint") +ax2.hlines( + y33, + xmin=0, + xmax=500, + color="purple", + linestyles="--", + linewidth=0.8, + label="Representative agent", +) +ax2.tick_params(axis="y", labelcolor=color2) plt.legend() # fig.suptitle('Zt and Market Price of Risk') fig.tight_layout() # otherwise the right y-label is slightly clipped -plt.savefig('Zt and interest rate' + '.png', dpi=500) +plt.savefig("Zt and interest rate" + ".png", dpi=500) plt.show() ####################################### @@ -110,22 +126,29 @@ y43 = f_parti_drop y44 = parti_drop fig, ax1 = plt.subplots(figsize=(15, 5)) -ax1.set_xlabel('Time in simulation, one random path') -ax1.set_ylabel('Zt', color=color5) +ax1.set_xlabel("Time in simulation, one random path") +ax1.set_ylabel("Zt", color=color5) ax1.plot(t, y0, color=color5, linewidth=0.5) -ax1.tick_params(axis='y', labelcolor=color5) +ax1.tick_params(axis="y", labelcolor=color5) ax2 = ax1.twinx() -ax2.set_ylabel('Consumption weighted bias', color=color2) +ax2.set_ylabel("Consumption weighted bias", color=color2) ax2.set_ylim([-0.5, 1]) -ax2.plot(t, y41, color='purple', linewidth=0.4, linestyle ='--', label='Market view, complete market') -ax2.plot(t, y42, color=color3, linewidth=0.4, label='Market view, no shorting') -ax2.plot(t, y43, color=color4, linewidth=0.4, label='Participant consumption share') -ax2.plot(t, y44, color=color2, linewidth=0.4, label='Participation rate') -ax2.tick_params(axis='y', labelcolor=color2) +ax2.plot( + t, + y41, + color="purple", + linewidth=0.4, + linestyle="--", + label="Market view, complete market", +) +ax2.plot(t, y42, color=color3, linewidth=0.4, label="Market view, no shorting") +ax2.plot(t, y43, color=color4, linewidth=0.4, label="Participant consumption share") +ax2.plot(t, y44, color=color2, linewidth=0.4, label="Participation rate") +ax2.tick_params(axis="y", labelcolor=color2) plt.legend() # fig.suptitle('Zt and market view') fig.tight_layout() # otherwise the right y-label is slightly clipped -plt.savefig('Zt and market bias' + '.png', dpi=500) +plt.savefig("Zt and market bias" + ".png", dpi=500) plt.show() # Regression @@ -138,66 +161,86 @@ invest_drop = pi_drop > 0 parti_rate_drop = invest_drop * cohort_size var = parti_rate_drop -y_label = 'participation composition of age groups' +y_label = "participation composition of age groups" y51 = np.sum(var[:, tau_cutoff1:], axis=1) y52 = np.sum(var[:, tau_cutoff2:], axis=1) y53 = np.sum(var[:, tau_cutoff3:], axis=1) -y54 = np.sum(var[:, ], axis=1) +y54 = np.sum( + var[ + :, + ], + axis=1, +) fig, ax1 = plt.subplots(figsize=(15, 5)) -ax1.set_ylabel('Participation composition of age groups', color=color2) +ax1.set_ylabel("Participation composition of age groups", color=color2) ax1.set_ylim([0, 1]) -ax1.fill_between(t, y51, color = 'steelblue', linewidth = 0.4, label = '20 < Age <= 35, youngest quartile') -ax1.fill_between(t, y52, y51, color = 'darkseagreen', linewidth = 0.4, label = '35 < Age <= 55') -ax1.fill_between(t, y53, y52, color = 'moccasin', linewidth = 0.4, label= '55 < Age <= 89') -ax1.fill_between(t, y54, y53, color = 'pink', linewidth = 0.4, label= 'Age > 89, oldest quartile') +ax1.fill_between( + t, y51, color="steelblue", linewidth=0.4, label="20 < Age <= 35, youngest quartile" +) +ax1.fill_between( + t, y52, y51, color="darkseagreen", linewidth=0.4, label="35 < Age <= 55" +) +ax1.fill_between(t, y53, y52, color="moccasin", linewidth=0.4, label="55 < Age <= 89") +ax1.fill_between( + t, y54, y53, color="pink", linewidth=0.4, label="Age > 89, oldest quartile" +) # ax1.fill_between(t, y51, color='lavender', linewidth=0.4, label='20 < Age <= 35, youngest quartile') # ax1.fill_between(t, y52, y51, color='lightsteelblue', linewidth=0.4, label='35 < Age <= 55') # ax1.fill_between(t, y53, y52, color='steelblue', linewidth=0.4, label='55 < Age <= 89') # ax1.fill_between(t, y54, y53, color='royalblue', linewidth=0.4, label='Age > 89, oldest quartile') -ax1.tick_params(axis='y', labelcolor=color2) +ax1.tick_params(axis="y", labelcolor=color2) plt.legend() ax2 = ax1.twinx() -ax2.set_xlabel('Time in simulation, one random path') -ax2.set_ylabel('Zt', color=color5) +ax2.set_xlabel("Time in simulation, one random path") +ax2.set_ylabel("Zt", color=color5) ax2.plot(t, y0, color=color5, linewidth=0.5) -ax2.tick_params(axis='y', labelcolor=color5) +ax2.tick_params(axis="y", labelcolor=color5) fig.tight_layout() # otherwise the right y-label is slightly clipped -plt.savefig('Zt and participation composition' + '.png', dpi = 500) +plt.savefig("Zt and participation composition" + ".png", dpi=500) plt.show() - # who want to short short_popu_drop = short_drop * cohort_size var = short_popu_drop y51 = np.sum(var[:, tau_cutoff1:], axis=1) y52 = np.sum(var[:, tau_cutoff2:], axis=1) y53 = np.sum(var[:, tau_cutoff3:], axis=1) -y54 = np.sum(var[:, ], axis=1) +y54 = np.sum( + var[ + :, + ], + axis=1, +) fig, ax1 = plt.subplots(figsize=(15, 5)) -ax1.set_ylabel('Short-sale demand composition of age groups', color=color2) +ax1.set_ylabel("Short-sale demand composition of age groups", color=color2) ax1.set_ylim([0, 0.02]) -ax1.fill_between(t, y51, color = 'steelblue', linewidth = 0.1, label = '20 < Age <= 35, youngest quartile') -ax1.fill_between(t, y52, y51, color = 'darkseagreen', linewidth = 0.1, label = '35 < Age <= 55') -ax1.fill_between(t, y53, y52, color = 'moccasin', linewidth = 0.1, label= '55 < Age <= 89') -ax1.fill_between(t, y54, y53, color = 'pink', linewidth = 0.1, label= 'Age > 89, oldest quartile') +ax1.fill_between( + t, y51, color="steelblue", linewidth=0.1, label="20 < Age <= 35, youngest quartile" +) +ax1.fill_between( + t, y52, y51, color="darkseagreen", linewidth=0.1, label="35 < Age <= 55" +) +ax1.fill_between(t, y53, y52, color="moccasin", linewidth=0.1, label="55 < Age <= 89") +ax1.fill_between( + t, y54, y53, color="pink", linewidth=0.1, label="Age > 89, oldest quartile" +) # ax1.fill_between(t, y51, color='lavender', linewidth=0.4, label='20 < Age <= 35, youngest quartile') # ax1.fill_between(t, y52, y51, color='lightsteelblue', linewidth=0.4, label='35 < Age <= 55') # ax1.fill_between(t, y53, y52, color='steelblue', linewidth=0.4, label='55 < Age <= 89') # ax1.fill_between(t, y54, y53, color='royalblue', linewidth=0.4, label='Age > 89, oldest quartile') -ax1.tick_params(axis='y', labelcolor=color2) +ax1.tick_params(axis="y", labelcolor=color2) plt.legend() ax2 = ax1.twinx() -ax2.set_xlabel('Time in simulation, one random path') -ax2.set_ylabel('Zt', color=color5) +ax2.set_xlabel("Time in simulation, one random path") +ax2.set_ylabel("Zt", color=color5) ax2.plot(t, y0, color=color5, linewidth=0.5) -ax2.tick_params(axis='y', labelcolor=color5) +ax2.tick_params(axis="y", labelcolor=color5) fig.tight_layout() # otherwise the right y-label is slightly clipped -#plt.savefig('Zt and short sale composition' + '.png', dpi = 500) +# plt.savefig('Zt and short sale composition' + '.png', dpi = 500) plt.show() - # illustrate who are investing and when do they quit nn = 10 length = len(t) @@ -205,7 +248,7 @@ starts = np.zeros(nn) for i in range(nn): start = int((i + 5) * 25 * (1 / dt)) - starts[i] = start*dt + starts[i] = start * dt for j in range(length): if j < start: pi_time_series[i, j] = np.nan @@ -214,30 +257,40 @@ a = pi_drop[j, cohort_rank] pi_time_series[i, j] = a if a == 0: - pi_time_series[i, j + 1: j + 8] = 0 - pi_time_series[i, j + 8:] = np.nan + pi_time_series[i, j + 1 : j + 8] = 0 + pi_time_series[i, j + 8 :] = np.nan break -colors = ['darkmagenta', 'midnightblue','green', 'saddlebrown', 'darkgreen','firebrick', 'purple', 'blue', 'olivedrab', 'darkviolet'] +colors = [ + "darkmagenta", + "midnightblue", + "green", + "saddlebrown", + "darkgreen", + "firebrick", + "purple", + "blue", + "olivedrab", + "darkviolet", +] fig, ax1 = plt.subplots(figsize=(15, 5)) -ax1.set_xlabel('Time in simulation, one random path') -ax1.set_ylabel('Zt', color=color5) +ax1.set_xlabel("Time in simulation, one random path") +ax1.set_ylabel("Zt", color=color5) ax1.plot(t, y0, color=color5, linewidth=0.5) -ax1.tick_params(axis='y', labelcolor=color5) +ax1.tick_params(axis="y", labelcolor=color5) ax2 = ax1.twinx() -ax2.set_ylabel('Investment in stock market', color=color2) +ax2.set_ylabel("Investment in stock market", color=color2) ax2.set_ylim([-0.1, 25]) for i in range(nn): y6 = pi_time_series[i] - plt.vlines(starts[i], ymax=25, ymin=0, color='grey', linestyle = '--', linewidth=0.4) + plt.vlines(starts[i], ymax=25, ymin=0, color="grey", linestyle="--", linewidth=0.4) ax2.plot(t, y6, color=colors[i], linewidth=0.4) -ax2.tick_params(axis='y', labelcolor=color2) +ax2.tick_params(axis="y", labelcolor=color2) fig.tight_layout() # otherwise the right y-label is slightly clipped -plt.savefig('Zt and pi time series' + '.png', dpi=500) +plt.savefig("Zt and pi time series" + ".png", dpi=500) plt.show() - ####################################### ############ GRAPH FIVE ############## ####################################### @@ -249,21 +302,21 @@ y71 = theta_free fig, ax1 = plt.subplots(figsize=(15, 5)) -ax1.set_xlabel('Time in simulation, one random path') -ax1.set_ylabel('Zt', color=color5) +ax1.set_xlabel("Time in simulation, one random path") +ax1.set_ylabel("Zt", color=color5) ax1.plot(t, y0, color=color5, linewidth=0.5) -ax1.tick_params(axis='y', labelcolor=color5) +ax1.tick_params(axis="y", labelcolor=color5) ax2 = ax1.twinx() -ax2.set_ylabel('Market price of risk', color=color2) +ax2.set_ylabel("Market price of risk", color=color2) ax2.set_ylim([-1, 1]) -ax2.plot(t, y21, color=color2, linewidth=0.4, label='Complete market') -ax2.plot(t, y22, color=color3, linewidth=0.4, label='Short-sale constraint') -ax2.plot(t, y71, color='magenta', linewidth=0.4, label='Rich can short') -#ax2.hlines(sigma_Y, xmin=0, xmax=500, color='purple', linestyles='--', linewidth=0.8, label='Representative agent') -ax2.tick_params(axis='y', labelcolor=color2) +ax2.plot(t, y21, color=color2, linewidth=0.4, label="Complete market") +ax2.plot(t, y22, color=color3, linewidth=0.4, label="Short-sale constraint") +ax2.plot(t, y71, color="magenta", linewidth=0.4, label="Rich can short") +# ax2.hlines(sigma_Y, xmin=0, xmax=500, color='purple', linestyles='--', linewidth=0.8, label='Representative agent') +ax2.tick_params(axis="y", labelcolor=color2) plt.legend() fig.tight_layout() # otherwise the right y-label is slightly clipped -plt.savefig('Zt and market price of risk, rich free' + '.png', dpi=500) +plt.savefig("Zt and market price of risk, rich free" + ".png", dpi=500) plt.show() @@ -274,110 +327,100 @@ y75 = Delta_bar_short_free fig, ax1 = plt.subplots(figsize=(15, 5)) -ax1.set_xlabel('Time in simulation, one random path') -ax1.set_ylabel('Zt', color=color5) +ax1.set_xlabel("Time in simulation, one random path") +ax1.set_ylabel("Zt", color=color5) ax1.plot(t, y0, color=color5, linewidth=0.5) -ax1.tick_params(axis='y', labelcolor=color5) +ax1.tick_params(axis="y", labelcolor=color5) ax2 = ax1.twinx() -ax2.set_ylabel('Market price of risk', color=color2) +ax2.set_ylabel("Market price of risk", color=color2) ax2.set_ylim([-0.5, 1]) -ax2.plot(t, y72, color='darkblue', linewidth=0.6, label='% investors') -ax2.plot(t, y73, color='darkgreen', linewidth=0.6, label='% Short sellers') -ax2.plot(t, y74, color='blue', linewidth=0.4, label='Average bias long') -ax2.plot(t, y75, color='magenta', linewidth=0.4, label='Average bias short') -#ax2.hlines(sigma_Y, xmin=0, xmax=500, color='purple', linestyles='--', linewidth=0.8, label='Representative agent') -ax2.tick_params(axis='y', labelcolor=color2) +ax2.plot(t, y72, color="darkblue", linewidth=0.6, label="% investors") +ax2.plot(t, y73, color="darkgreen", linewidth=0.6, label="% Short sellers") +ax2.plot(t, y74, color="blue", linewidth=0.4, label="Average bias long") +ax2.plot(t, y75, color="magenta", linewidth=0.4, label="Average bias short") +# ax2.hlines(sigma_Y, xmin=0, xmax=500, color='purple', linestyles='--', linewidth=0.8, label='Representative agent') +ax2.tick_params(axis="y", labelcolor=color2) plt.legend() fig.tight_layout() # otherwise the right y-label is slightly clipped -plt.savefig('Zt and participation, rich free' + '.png', dpi=500) +plt.savefig("Zt and participation, rich free" + ".png", dpi=500) plt.show() - ####################################### ############ GRAPH SIX ############### ####################################### # test if subjective risk premia comove less with shocks / cyclicality of perceived risk premia -theta_hat = theta_drop + Delta_bar_parti_drop # wealth weighted average subjective risk premia +theta_hat = ( + theta_drop + Delta_bar_parti_drop +) # wealth weighted average subjective risk premia # regression: -x1 = dZ.reshape((-1,1)) +x1 = dZ.reshape((-1, 1)) y_obj_theta = theta_drop y_subj_theta = theta_hat model_obj_theta = LinearRegression().fit(x1, y_obj_theta) model_subj_theta = LinearRegression().fit(x1, y_subj_theta) - - - - - - - - - - - y2 = popu_long_collect y22 = popu_long_free y23 = theta_comp y3 = age y4 = n_parti -tail = '_rich_free' +tail = "_rich_free" fig, ax1 = plt.subplots(figsize=(10, 5)) -ax1.set_xlabel('Time in simulation, one random path') -ax1.set_ylabel('Zt', color=color5) +ax1.set_xlabel("Time in simulation, one random path") +ax1.set_ylabel("Zt", color=color5) ax1.plot(t, y1, color=color5, linewidth=0.5) -ax1.tick_params(axis='y', labelcolor=color5) +ax1.tick_params(axis="y", labelcolor=color5) ax2 = ax1.twinx() # ax2.set_ylabel('Population investing in stocks', color = color2) -ax2.set_ylabel('Bias in belief and learning', color=color2) +ax2.set_ylabel("Bias in belief and learning", color=color2) ax2.set_ylim([-0.5, 0.5]) ax2.plot(t, y22, color=color2, linewidth=0.4) ax2.plot(t, y23, color=color3, linewidth=0.4) ax2.plot(t, y24, color=color4, linewidth=0.4) -ax2.tick_params(axis='y', labelcolor=color2) +ax2.tick_params(axis="y", labelcolor=color2) # fig.suptitle('Zt and Shorting Rate') fig.tight_layout() # otherwise the right y-label is slightly clipped # plt.savefig('Zt and Shorting Rate Comparison' + tail + '.jpg') -plt.savefig('Zt and bias time series' + '.png', dpi=500) +plt.savefig("Zt and bias time series" + ".png", dpi=500) plt.show() fig, ax1 = plt.subplots() -color1 = 'r' -ax1.set_xlabel('time in simulation, one random path') -ax1.set_ylabel('Zt', color=color1) +color1 = "r" +ax1.set_xlabel("time in simulation, one random path") +ax1.set_ylabel("Zt", color=color1) ax1.plot(t, y1, color=color1, linewidth=0.8) -ax1.tick_params(axis='y', labelcolor=color1) +ax1.tick_params(axis="y", labelcolor=color1) ax2 = ax1.twinx() -color2 = 'b' -ax2.set_ylabel('Average age holding stocks', color=color2) +color2 = "b" +ax2.set_ylabel("Average age holding stocks", color=color2) ax2.set_ylim([0, 100]) ax2.plot(t, y3, color=color2, linewidth=0.8) -ax2.tick_params(axis='y', labelcolor=color2) -fig.suptitle('Zt and Average Participant Age') +ax2.tick_params(axis="y", labelcolor=color2) +fig.suptitle("Zt and Average Participant Age") fig.tight_layout() # otherwise the right y-label is slightly clipped -plt.savefig('Zt and Average Participant Age' + tail + '.jpg') +plt.savefig("Zt and Average Participant Age" + tail + ".jpg") plt.show() fig, ax1 = plt.subplots() -color1 = 'r' -ax1.set_xlabel('time in simulation, one random path') -ax1.set_ylabel('Zt', color=color1) +color1 = "r" +ax1.set_xlabel("time in simulation, one random path") +ax1.set_ylabel("Zt", color=color1) ax1.plot(t, y1, color=color1, linewidth=0.8) -ax1.tick_params(axis='y', labelcolor=color1) +ax1.tick_params(axis="y", labelcolor=color1) ax2 = ax1.twinx() -color2 = 'b' -ax2.set_ylabel('% of cohorts holding stocks', color=color2) +color2 = "b" +ax2.set_ylabel("% of cohorts holding stocks", color=color2) ax2.set_ylim([0, 1]) ax2.plot(t, y4, color=color2, linewidth=0.8) -ax2.tick_params(axis='y', labelcolor=color2) -fig.suptitle('Zt and % of Cohorts Participate') +ax2.tick_params(axis="y", labelcolor=color2) +fig.suptitle("Zt and % of Cohorts Participate") fig.tight_layout() # otherwise the right y-label is slightly clipped -plt.savefig('Zt and Cohorts Participate' + tail + '.jpg') +plt.savefig("Zt and Cohorts Participate" + tail + ".jpg") plt.show() ############################## @@ -387,54 +430,54 @@ y4 = n_parti_drop fig, ax1 = plt.subplots() -color1 = 'r' -ax1.set_xlabel('time in simulation, one random path') -ax1.set_ylabel('Zt', color=color1) +color1 = "r" +ax1.set_xlabel("time in simulation, one random path") +ax1.set_ylabel("Zt", color=color1) ax1.plot(t, y1, color=color1, linewidth=0.8) -ax1.tick_params(axis='y', labelcolor=color1) +ax1.tick_params(axis="y", labelcolor=color1) ax2 = ax1.twinx() -color2 = 'b' -ax2.set_ylabel('Population holding stocks, drop', color=color2) +color2 = "b" +ax2.set_ylabel("Population holding stocks, drop", color=color2) ax2.set_ylim([0, 1]) ax2.plot(t, y2, color=color2, linewidth=0.8) -ax2.tick_params(axis='y', labelcolor=color2) -fig.suptitle('Zt and Participation Rate, if Drop') +ax2.tick_params(axis="y", labelcolor=color2) +fig.suptitle("Zt and Participation Rate, if Drop") fig.tight_layout() # otherwise the right y-label is slightly clipped -plt.savefig('Zt and Participation Rate' + '_drop' + '.jpg') +plt.savefig("Zt and Participation Rate" + "_drop" + ".jpg") plt.show() fig, ax1 = plt.subplots() -color1 = 'r' -ax1.set_xlabel('time in simulation, one random path') -ax1.set_ylabel('Zt', color=color1) +color1 = "r" +ax1.set_xlabel("time in simulation, one random path") +ax1.set_ylabel("Zt", color=color1) ax1.plot(t, y1, color=color1, linewidth=0.8) -ax1.tick_params(axis='y', labelcolor=color1) +ax1.tick_params(axis="y", labelcolor=color1) ax2 = ax1.twinx() -color2 = 'b' -ax2.set_ylabel('Average age holding stocks', color=color2) +color2 = "b" +ax2.set_ylabel("Average age holding stocks", color=color2) ax2.set_ylim([0, 100]) ax2.plot(t, y3, color=color2, linewidth=0.8) -ax2.tick_params(axis='y', labelcolor=color2) -fig.suptitle('Zt and Average Participant Age, if Drop') +ax2.tick_params(axis="y", labelcolor=color2) +fig.suptitle("Zt and Average Participant Age, if Drop") fig.tight_layout() # otherwise the right y-label is slightly clipped -plt.savefig('Zt and Average Participant Age' + '_drop' + '.jpg') +plt.savefig("Zt and Average Participant Age" + "_drop" + ".jpg") plt.show() fig, ax1 = plt.subplots() -color1 = 'r' -ax1.set_xlabel('time in simulation, one random path') -ax1.set_ylabel('Zt', color=color1) +color1 = "r" +ax1.set_xlabel("time in simulation, one random path") +ax1.set_ylabel("Zt", color=color1) ax1.plot(t, y1, color=color1, linewidth=0.8) -ax1.tick_params(axis='y', labelcolor=color1) +ax1.tick_params(axis="y", labelcolor=color1) ax2 = ax1.twinx() -color2 = 'b' -ax2.set_ylabel('% of cohorts holding stocks', color=color2) +color2 = "b" +ax2.set_ylabel("% of cohorts holding stocks", color=color2) ax2.set_ylim([0, 1]) ax2.plot(t, y4, color=color2, linewidth=0.8) -ax2.tick_params(axis='y', labelcolor=color2) -fig.suptitle('Zt and % of Cohorts Participate, if Drop') +ax2.tick_params(axis="y", labelcolor=color2) +fig.suptitle("Zt and % of Cohorts Participate, if Drop") fig.tight_layout() # otherwise the right y-label is slightly clipped -plt.savefig('Zt and Cohorts Participate' + '_drop' + '.jpg') +plt.savefig("Zt and Cohorts Participate" + "_drop" + ".jpg") plt.show() y1 = Z @@ -443,54 +486,54 @@ y4 = n_parti_comp fig, ax1 = plt.subplots() -color1 = 'r' -ax1.set_xlabel('time in simulation, one random path') -ax1.set_ylabel('Zt', color=color1) +color1 = "r" +ax1.set_xlabel("time in simulation, one random path") +ax1.set_ylabel("Zt", color=color1) ax1.plot(t, y1, color=color1, linewidth=0.8) -ax1.tick_params(axis='y', labelcolor=color1) +ax1.tick_params(axis="y", labelcolor=color1) ax2 = ax1.twinx() -color2 = 'b' -ax2.set_ylabel('Population longing stocks', color=color2) +color2 = "b" +ax2.set_ylabel("Population longing stocks", color=color2) ax2.set_ylim([0, 1]) ax2.plot(t, y2, color=color2, linewidth=0.8) -ax2.tick_params(axis='y', labelcolor=color2) -fig.suptitle('Zt and Participation Rate, complete market') +ax2.tick_params(axis="y", labelcolor=color2) +fig.suptitle("Zt and Participation Rate, complete market") fig.tight_layout() # otherwise the right y-label is slightly clipped -plt.savefig('Zt and Participation Rate' + '_comp' + '.jpg') +plt.savefig("Zt and Participation Rate" + "_comp" + ".jpg") plt.show() fig, ax1 = plt.subplots() -color1 = 'r' -ax1.set_xlabel('time in simulation, one random path') -ax1.set_ylabel('Zt', color=color1) +color1 = "r" +ax1.set_xlabel("time in simulation, one random path") +ax1.set_ylabel("Zt", color=color1) ax1.plot(t, y1, color=color1, linewidth=0.8) -ax1.tick_params(axis='y', labelcolor=color1) +ax1.tick_params(axis="y", labelcolor=color1) ax2 = ax1.twinx() -color2 = 'b' -ax2.set_ylabel('Average age longing stocks', color=color2) +color2 = "b" +ax2.set_ylabel("Average age longing stocks", color=color2) ax2.set_ylim([0, 100]) ax2.plot(t, y3, color=color2, linewidth=0.8) -ax2.tick_params(axis='y', labelcolor=color2) -fig.suptitle('Zt and Average Participant Age, complete market') +ax2.tick_params(axis="y", labelcolor=color2) +fig.suptitle("Zt and Average Participant Age, complete market") fig.tight_layout() # otherwise the right y-label is slightly clipped -plt.savefig('Zt and Average Participant Age' + '_comp' + '.jpg') +plt.savefig("Zt and Average Participant Age" + "_comp" + ".jpg") plt.show() fig, ax1 = plt.subplots() -color1 = 'r' -ax1.set_xlabel('time in simulation, one random path') -ax1.set_ylabel('Zt', color=color1) +color1 = "r" +ax1.set_xlabel("time in simulation, one random path") +ax1.set_ylabel("Zt", color=color1) ax1.plot(t, y1, color=color1, linewidth=0.8) -ax1.tick_params(axis='y', labelcolor=color1) +ax1.tick_params(axis="y", labelcolor=color1) ax2 = ax1.twinx() -color2 = 'b' -ax2.set_ylabel('% of cohorts longing stocks', color=color2) +color2 = "b" +ax2.set_ylabel("% of cohorts longing stocks", color=color2) ax2.set_ylim([0, 1]) ax2.plot(t, y4, color=color2, linewidth=0.8) -ax2.tick_params(axis='y', labelcolor=color2) -fig.suptitle('Zt and % of Cohorts Participate, complete market') +ax2.tick_params(axis="y", labelcolor=color2) +fig.suptitle("Zt and % of Cohorts Participate, complete market") fig.tight_layout() # otherwise the right y-label is slightly clipped -plt.savefig('Zt and Cohorts Participate' + '_comp' + '.jpg') +plt.savefig("Zt and Cohorts Participate" + "_comp" + ".jpg") plt.show() ############################# @@ -501,22 +544,22 @@ y3 = r_comp fig, ax1 = plt.subplots() -color1 = 'r' -ax1.set_xlabel('time in simulation, one random path') -ax1.set_ylabel('Zt', color=color1) +color1 = "r" +ax1.set_xlabel("time in simulation, one random path") +ax1.set_ylabel("Zt", color=color1) ax1.plot(t, y1, color=color1, linewidth=0.8) -ax1.tick_params(axis='y', labelcolor=color1) +ax1.tick_params(axis="y", labelcolor=color1) ax2 = ax1.twinx() -color2 = 'b' -color3 = 'g' -ax2.set_ylabel('Interest Rate', color=color2) +color2 = "b" +color3 = "g" +ax2.set_ylabel("Interest Rate", color=color2) ax2.set_ylim([0, 0.05]) ax2.plot(t, y2, color=color2, linewidth=0.5) ax2.plot(t, y3, color=color3, linewidth=0.5) -ax2.tick_params(axis='y', labelcolor=color2) -fig.suptitle('Zt and Interest Rate') +ax2.tick_params(axis="y", labelcolor=color2) +fig.suptitle("Zt and Interest Rate") fig.tight_layout() # otherwise the right y-label is slightly clipped -plt.savefig('Zt and interest rate' + '.jpg') +plt.savefig("Zt and interest rate" + ".jpg") plt.show() y2 = theta_drop @@ -524,24 +567,24 @@ y23 = theta_comp fig, ax1 = plt.subplots() -color1 = 'r' -ax1.set_xlabel('time in simulation, one random path') -ax1.set_ylabel('Zt', color=color1) +color1 = "r" +ax1.set_xlabel("time in simulation, one random path") +ax1.set_ylabel("Zt", color=color1) ax1.plot(t, y1, color=color1, linewidth=0.8) -ax1.tick_params(axis='y', labelcolor=color1) +ax1.tick_params(axis="y", labelcolor=color1) ax2 = ax1.twinx() -color2 = 'b' -color3 = 'g' -color4 = 'm' -ax2.set_ylabel('equity risk premium', color=color2) +color2 = "b" +color3 = "g" +color4 = "m" +ax2.set_ylabel("equity risk premium", color=color2) ax2.set_ylim([-0.5, 1]) ax2.plot(t, y2, color=color2, linewidth=0.5) ax2.plot(t, y22, color=color3, linewidth=0.5) ax2.plot(t, y23, color=color4, linewidth=0.5) -ax2.tick_params(axis='y', labelcolor=color2) -fig.suptitle('Zt and Equity Risk Premium, drop=b, free=g, comp=m') +ax2.tick_params(axis="y", labelcolor=color2) +fig.suptitle("Zt and Equity Risk Premium, drop=b, free=g, comp=m") fig.tight_layout() # otherwise the right y-label is slightly clipped -plt.savefig('Zt and equity risk premium' + '.jpg') +plt.savefig("Zt and equity risk premium" + ".jpg") plt.show() # f_dt_drop = f_drop * dt @@ -557,7 +600,7 @@ parti_rate_drop = invest_drop * cohort_size var = f_dt_free -y_label = 'wealth share of age groups' +y_label = "wealth share of age groups" # y_label = 'bias of age groups' # y_label = 'participation rate of age groups' # y_label = 'participation composition of age groups' @@ -568,24 +611,24 @@ y5 = np.sum(var[:, :tau_cutoff3], axis=1) fig, ax1 = plt.subplots() -color1 = 'r' -ax1.set_xlabel('time in simulation, one random path') -ax1.set_ylabel('Zt', color=color1) +color1 = "r" +ax1.set_xlabel("time in simulation, one random path") +ax1.set_ylabel("Zt", color=color1) ax1.plot(t, y1, color=color1, linewidth=0.8) -ax1.tick_params(axis='y', labelcolor=color1) +ax1.tick_params(axis="y", labelcolor=color1) ax2 = ax1.twinx() -color2 = 'b' -color3 = 'g' -color4 = 'y' -color5 = 'm' +color2 = "b" +color3 = "g" +color4 = "y" +color5 = "m" ax2.set_ylabel(y_label, color=color2) ax2.set_ylim([0, 1]) ax2.plot(t, y2, color=color2, linewidth=0.5) ax2.plot(t, y3, color=color3, linewidth=0.5) ax2.plot(t, y4, color=color4, linewidth=0.5) ax2.plot(t, y5, color=color5, linewidth=0.5) -ax2.tick_params(axis='y', labelcolor=color2) -fig.suptitle('Zt and ' + y_label + '(blue, green, yellow, magenta)') +ax2.tick_params(axis="y", labelcolor=color2) +fig.suptitle("Zt and " + y_label + "(blue, green, yellow, magenta)") fig.tight_layout() # otherwise the right y-label is slightly clipped -plt.savefig('Zt and ' + y_label + tail + '.jpg') +plt.savefig("Zt and " + y_label + tail + ".jpg") plt.show() diff --git a/main.py b/main.py index ad22803..a750337 100644 --- a/main.py +++ b/main.py @@ -16,10 +16,12 @@ # modes = ['keep', 'drop', 'comp', 'rich_free', 'back_collect', 'back_renew'] # modes = ['drop'] # modes = ['rich_free'] -modes = ['keep'] +modes = ["keep"] # The main loop builds up the economy with a large number of cohorts, and simulates the stationary economy forward -survey_view_parti_matrix = np.zeros((Mpaths, Nt)) # average perceived risk premia among investors +survey_view_parti_matrix = np.zeros( + (Mpaths, Nt) +) # average perceived risk premia among investors survey_view_long_only_matrix = np.zeros((Mpaths, Nt)) survey_view_can_short_matrix = np.zeros((Mpaths, Nt)) survey_view_parti_young_matrix = np.zeros((Mpaths, Nt)) @@ -39,14 +41,18 @@ time_s = time.time() if k == Mpaths - 1: # in the last round, use the shocks seen in the slides - dZ_build = np.load('dZt_build_demo.npy') # dZt for the build function - dZ = np.load('dZt_demo.npy') # dZt for the simulate function + dZ_build = np.load("dZt_build_demo.npy") # dZt for the build function + dZ = np.load("dZt_demo.npy") # dZt for the simulate function else: - dZ_build = dt ** 0.5 * np.random.randn(int(Nc - 1)) # dZt for the build function - dZ = dt ** 0.5 * np.random.randn(Nt) # dZt for the simulate function + dZ_build = dt**0.5 * np.random.randn( + int(Nc - 1) + ) # dZt for the build function + dZ = dt**0.5 * np.random.randn(Nt) # dZt for the simulate function - biasvec = dZ_build[-Npre:] # dZt in the building cohorts stage, but also used in the simulation function + biasvec = dZ_build[ + -Npre: + ] # dZt in the building cohorts stage, but also used in the simulation function dY_Y_matrix[k, :] = mu_Y * dt + sigma_Y * dZ @@ -55,7 +61,7 @@ # Z_matrix[k, :] = Z for mode in modes: - if mode == 'keep' or mode == 'drop' or mode == 'complete': + if mode == "keep" or mode == "drop" or mode == "complete": ( mu_S, mu_S_s, @@ -73,14 +79,30 @@ w_cohort, age_parti, n_parti, - ) = simulate(mode, Nc, Nt, dt, rho, nu, Vhat, mu_Y, sigma_Y, tax, beta, Npre, T_hat, dZ_build, dZ, tau, - cohort_size) + ) = simulate( + mode, + Nc, + Nt, + dt, + rho, + nu, + Vhat, + mu_Y, + sigma_Y, + tax, + beta, + Npre, + T_hat, + dZ_build, + dZ, + tau, + cohort_size, + ) invest_tracker = pi > 0 theta_mat = np.transpose(np.tile(theta, (Nc, 1))) - - if mode == 'rich_free' or mode == 'back_collect' or mode == 'back_renew': + if mode == "rich_free" or mode == "back_collect" or mode == "back_renew": ( mu_S, mu_S_s, @@ -111,17 +133,36 @@ Delta_bar_parti, Delta_bar_long, Delta_bar_short, - ) = simulate_partial_constraint(mode, Nc, Nt, dt, rho, nu, Vhat, mu_Y, sigma_Y, tax, beta, Npre, - T_hat, dZ_build, dZ, tau, cohort_size) + ) = simulate_partial_constraint( + mode, + Nc, + Nt, + dt, + rho, + nu, + Vhat, + mu_Y, + sigma_Y, + tax, + beta, + Npre, + T_hat, + dZ_build, + dZ, + tau, + cohort_size, + ) theta_mat = np.transpose(np.tile(theta, (Nc, 1))) long_only_weights = (invest_tracker - can_short_tracker) * cohort_size can_short_weights = can_short_tracker * cohort_size - survey_view_long_only_matrix[k] = np.average((Delta + theta_mat) * sigma_S, - weights=long_only_weights, axis=1) - survey_view_can_short_matrix[k] = np.average((Delta + theta_mat) * sigma_S, - weights=can_short_weights, axis=1) + survey_view_long_only_matrix[k] = np.average( + (Delta + theta_mat) * sigma_S, weights=long_only_weights, axis=1 + ) + survey_view_can_short_matrix[k] = np.average( + (Delta + theta_mat) * sigma_S, weights=can_short_weights, axis=1 + ) erp_S = mu_S - r # erp_S_s = mu_S_s - np.reshape(r, (Nt, 1)) @@ -164,28 +205,39 @@ # todo: correct the code about survey view cohort_size_mat = np.tile(cohort_size, (Nt, 1)) survey_view_parti = (Delta + theta_mat) * invest_tracker * sigma_S - survey_view_parti_matrix[k] = np.average(survey_view_parti, - weights=cohort_size_mat, axis=1) + survey_view_parti_matrix[k] = np.average( + survey_view_parti, weights=cohort_size_mat, axis=1 + ) - weights_zero = (np.sum(survey_view_parti[:, tau_cutoff1:], axis=1) == 0) + weights_zero = np.sum(survey_view_parti[:, tau_cutoff1:], axis=1) == 0 view_copy = np.copy(survey_view_parti) a = np.where(weights_zero == 1) view_copy[a, :] = np.nan - survey_view_parti_young_matrix[k] = np.average(view_copy[:, tau_cutoff1:], - weights = cohort_size_mat[:, tau_cutoff1:], axis=1) + survey_view_parti_young_matrix[k] = np.average( + view_copy[:, tau_cutoff1:], weights=cohort_size_mat[:, tau_cutoff1:], axis=1 + ) - weights_zero1 = (np.sum(survey_view_parti[:, tau_cutoff3:tau_cutoff2], axis=1) == 0) + weights_zero1 = ( + np.sum(survey_view_parti[:, tau_cutoff3:tau_cutoff2], axis=1) == 0 + ) view_copy1 = np.copy(survey_view_parti) a1 = np.where(weights_zero1 == 1) view_copy1[a1, :] = np.nan - survey_view_parti_old_matrix[k] = np.average(view_copy1[:, tau_cutoff3:tau_cutoff2], - weights=cohort_size_mat[:, tau_cutoff3:tau_cutoff2], axis=1) + survey_view_parti_old_matrix[k] = np.average( + view_copy1[:, tau_cutoff3:tau_cutoff2], + weights=cohort_size_mat[:, tau_cutoff3:tau_cutoff2], + axis=1, + ) obj_rp_matrix[k] = theta * sigma_S parti_track = cohort_size_mat * invest_tracker - popu_parti_young_matrix[k] = np.sum(parti_track[:, tau_cutoff1:], axis=1) # the first age quartile - popu_parti_old_matrix[k] = np.sum(parti_track[:, tau_cutoff3:tau_cutoff2], axis=1) # the third age quartile + popu_parti_young_matrix[k] = np.sum( + parti_track[:, tau_cutoff1:], axis=1 + ) # the first age quartile + popu_parti_old_matrix[k] = np.sum( + parti_track[:, tau_cutoff3:tau_cutoff2], axis=1 + ) # the third age quartile # # for graphs specific to one random path: # if dZ_build == np.load('dZt_build_demo.npy'): @@ -471,9 +523,9 @@ for j, horizon in enumerate(horizons): if i == 0: - header_j = str(horizon) + '-month' + header_j = str(horizon) + "-month" header.append(header_j) - x = (x_path[horizon:-horizon] - x_path[:-horizon * 2]) / (horizon * dt) + x = (x_path[horizon:-horizon] - x_path[: -horizon * 2]) / (horizon * dt) x = x / np.std(x) x = x.reshape(-1, 1) x = sm.add_constant(x) @@ -503,21 +555,23 @@ result_sub = np.mean(results_sub_matrix, axis=0) # to table: -index = ['coef', 't-stats', 'R2'] +index = ["coef", "t-stats", "R2"] n = len(index) for i in range(2): reg_data = np.empty((n, m)) var = result_obj if i == 0 else result_sub for j in range(n): - reg_data[j] = var[:,j] - print('result_obj' if i == 0 else 'result_sub') - print(tabulate.tabulate(reg_data, headers=header, showindex=index, floatfmt=".4f", tablefmt='fancy_grid')) - - - - - - + reg_data[j] = var[:, j] + print("result_obj" if i == 0 else "result_sub") + print( + tabulate.tabulate( + reg_data, + headers=header, + showindex=index, + floatfmt=".4f", + tablefmt="fancy_grid", + ) + ) ####################################### @@ -567,7 +621,12 @@ start_t = 0 -x_variables = [popu_parti_matrix, popu_parti_young_matrix, popu_parti_old_matrix, age_parti_matrix] +x_variables = [ + popu_parti_matrix, + popu_parti_young_matrix, + popu_parti_old_matrix, + age_parti_matrix, +] m = len(x_variables) # np.cumsum(dR_matrix[i]) coeff_matrix1 = np.zeros((Mpaths, m, 3)) @@ -620,9 +679,18 @@ reg_pvalues1 = np.average(pvalue_matrix1, axis=0) reg_tstats1 = np.average(tstats_matrix1, axis=0) reg_rsqrd1 = np.average(rsqrd_matrix1, axis=0) -reg_data = np.empty((5,8)) -header = ['(1) parti rate', '(2)', '(3) parti rate, young', '(4)', '(5) parti rate, old', '(6)', '(7) parti age', '(8)'] -index = ['coef', 't-stats', 'coef_x2', 't-stats_x2', 'R2'] +reg_data = np.empty((5, 8)) +header = [ + "(1) parti rate", + "(2)", + "(3) parti rate, young", + "(4)", + "(5) parti rate, old", + "(6)", + "(7) parti age", + "(8)", +] +index = ["coef", "t-stats", "coef_x2", "t-stats_x2", "R2"] for i in range(4): reg_data[0, i * 2] = reg_coeffs1[i, 0] reg_data[1, i * 2] = reg_tstats1[i, 0] @@ -636,13 +704,21 @@ reg_data[3, i * 2 + 1] = reg_tstats1[i, 2] reg_data[4, i * 2 + 1] = reg_rsqrd1[i, 1] -print(tabulate.tabulate(reg_data, headers=header, showindex = index, floatfmt=".4f", tablefmt='fancy_grid')) - +print( + tabulate.tabulate( + reg_data, headers=header, showindex=index, floatfmt=".4f", tablefmt="fancy_grid" + ) +) #### horizons = [1, 3, 6, 12, 36, 60, 120] -x_variables = [popu_parti_matrix, popu_parti_young_matrix, popu_parti_old_matrix, age_parti_matrix] +x_variables = [ + popu_parti_matrix, + popu_parti_young_matrix, + popu_parti_old_matrix, + age_parti_matrix, +] # x_variables = [popu_parti_matrix, popu_parti_young_matrix, popu_parti_old_matrix, age_parti_matrix] m = len(x_variables) n = len(horizons) @@ -657,14 +733,16 @@ path_r = np.cumsum(r_matrix[i]) path_x2 = survey_view_parti_matrix[i] path_x2 = path_x2 / np.std(path_x2) - for j,horizon in enumerate(horizons): - y_raw = (path_y[start_t + horizon + 1:] - path_y[start_t + 1: -horizon]) / (horizon * dt) - #dR is the return from t-1 to t, and thus have to move 1 to have returns from t to t+1 + for j, horizon in enumerate(horizons): + y_raw = (path_y[start_t + horizon + 1 :] - path_y[start_t + 1 : -horizon]) / ( + horizon * dt + ) + # dR is the return from t-1 to t, and thus have to move 1 to have returns from t to t+1 y_predict = y_raw.reshape(-1, 1) for k, var in enumerate(x_variables): path_x = var[i] # univariate regressions - x1_raw = path_x[start_t: -1-horizon] # participation rate at time t + x1_raw = path_x[start_t : -1 - horizon] # participation rate at time t x1_raw = x1_raw / np.std(x1_raw) x1_lag = x1_raw.reshape(-1, 1) x1_lag2 = sm.add_constant(x1_lag) @@ -677,7 +755,9 @@ rsqrd_matrix2[i, j, k, 0] = est.rsquared # bivariate regressions - x2_lag = path_x2[start_t: -1-horizon] # average perceived risk premia at time t + x2_lag = path_x2[ + start_t : -1 - horizon + ] # average perceived risk premia at time t x_lag = np.append(x1_lag, x2_lag) x_lag = np.transpose(x_lag.reshape(2, -1)) @@ -703,9 +783,17 @@ for j in range(n): horizon = horizons[j] reg_data = np.empty((5, 8)) - header = ['(1) parti rate', '(2)', '(3) parti rate, young', '(4)', '(5) parti rate, old', '(6)', '(7) parti age', - '(8)'] - index = ['coef', 't-stats', 'coef_x2', 't-stats_x2', 'R2'] + header = [ + "(1) parti rate", + "(2)", + "(3) parti rate, young", + "(4)", + "(5) parti rate, old", + "(6)", + "(7) parti age", + "(8)", + ] + index = ["coef", "t-stats", "coef_x2", "t-stats_x2", "R2"] for i in range(4): reg_data[0, i * 2] = reg_coeffs2[j, i, 0] reg_data[1, i * 2] = reg_tstats2[j, i, 0] @@ -719,15 +807,26 @@ reg_data[3, i * 2 + 1] = reg_tstats2[j, i, 2] reg_data[4, i * 2 + 1] = reg_rsqrd2[j, i, 1] - print(str(horizon) + '-month') - print(tabulate.tabulate(reg_data, headers=header, showindex=index, floatfmt=".4f", tablefmt='fancy_grid')) - - + print(str(horizon) + "-month") + print( + tabulate.tabulate( + reg_data, + headers=header, + showindex=index, + floatfmt=".4f", + tablefmt="fancy_grid", + ) + ) #### horizons = [1, 3, 6, 12, 36, 60, 120] -x_variables = [popu_parti_matrix, popu_parti_young_matrix, popu_parti_old_matrix, age_parti_matrix] +x_variables = [ + popu_parti_matrix, + popu_parti_young_matrix, + popu_parti_old_matrix, + age_parti_matrix, +] # x_variables = [popu_parti_matrix, popu_parti_young_matrix, popu_parti_old_matrix, age_parti_matrix] m = len(x_variables) n = len(horizons) @@ -742,15 +841,17 @@ path_r = np.cumsum(r_matrix[i]) path_x2 = survey_view_parti_matrix[i] path_x2 = path_x2 / np.std(path_x2) - for j,horizon in enumerate(horizons): + for j, horizon in enumerate(horizons): # y_raw = path_y[start_t + horizon + 1:] - path_y[start_t + 1: -horizon] #dR is the return from t-1 to t, and thus have to move 1 to have returns from t to t+1 - y_raw = ((path_y[start_t + horizon + 1:] - path_y[start_t + 1: -horizon]) \ - - (path_r[start_t + horizon: -1] - path_r[start_t : -horizon - 1])) / (horizon * dt) + y_raw = ( + (path_y[start_t + horizon + 1 :] - path_y[start_t + 1 : -horizon]) + - (path_r[start_t + horizon : -1] - path_r[start_t : -horizon - 1]) + ) / (horizon * dt) y_predict = y_raw.reshape(-1, 1) for k, var in enumerate(x_variables): path_x = var[i] # univariate regressions - x1_raw = path_x[start_t: -1-horizon] # participation rate at time t + x1_raw = path_x[start_t : -1 - horizon] # participation rate at time t x1_raw = x1_raw / np.std(x1_raw) x1_lag = x1_raw.reshape(-1, 1) x1_lag2 = sm.add_constant(x1_lag) @@ -763,7 +864,7 @@ rsqrd_matrix2[i, j, k, 0] = est.rsquared # bivariate regressions - x2_lag = path_x2[start_t: -horizon - 1] + x2_lag = path_x2[start_t : -horizon - 1] x_lag = np.append(x1_lag, x2_lag) x_lag = np.transpose(x_lag.reshape(2, -1)) @@ -789,9 +890,17 @@ for j in range(n): horizon = horizons[j] reg_data = np.empty((5, 8)) - header = ['(1) parti rate', '(2)', '(3) parti rate, young', '(4)', '(5) parti rate, old', '(6)', '(7) parti age', - '(8)'] - index = ['coef', 't-stats', 'coef_x2', 't-stats_x2', 'R2'] + header = [ + "(1) parti rate", + "(2)", + "(3) parti rate, young", + "(4)", + "(5) parti rate, old", + "(6)", + "(7) parti age", + "(8)", + ] + index = ["coef", "t-stats", "coef_x2", "t-stats_x2", "R2"] for i in range(4): reg_data[0, i * 2] = reg_coeffs2[j, i, 0] reg_data[1, i * 2] = reg_tstats2[j, i, 0] @@ -805,5 +914,13 @@ reg_data[3, i * 2 + 1] = reg_tstats2[j, i, 2] reg_data[4, i * 2 + 1] = reg_rsqrd2[j, i, 1] - print(str(horizon) + '-month') - print(tabulate.tabulate(reg_data, headers=header, showindex=index, floatfmt=".4f", tablefmt='fancy_grid')) \ No newline at end of file + print(str(horizon) + "-month") + print( + tabulate.tabulate( + reg_data, + headers=header, + showindex=index, + floatfmt=".4f", + tablefmt="fancy_grid", + ) + ) diff --git a/src/RES tax.py b/src/RES tax.py index 473cad2..cdbbdd2 100644 --- a/src/RES tax.py +++ b/src/RES tax.py @@ -7,9 +7,9 @@ ################# Define functions #################################################################### def PostVar(sigY, Vbar, tau): if type(tau) == np.ndarray: - V = sigY ** 2 * Vbar / (sigY ** 2 * np.ones(len(tau)) + Vbar * tau) + V = sigY**2 * Vbar / (sigY**2 * np.ones(len(tau)) + Vbar * tau) else: - V = sigY ** 2 * Vbar / (sigY ** 2 + Vbar * tau) + V = sigY**2 * Vbar / (sigY**2 + Vbar * tau) return V @@ -17,7 +17,7 @@ def BuildUpCohortsMAIN(dZt, Nt, dt, rho, nu, Vbar, mu_Y, sigma_Y, bet, That): # builds up a sufficiently large set of cohorts Npre = int(That / dt) # Number of pre-trading observations Zt = np.insert(np.cumsum(dZt), 0, 0) - yg = (mu_Y - 0.5 * sigma_Y ** 2) * dt * np.ones(int(Nt - 1)) + sigma_Y * dZt + yg = (mu_Y - 0.5 * sigma_Y**2) * dt * np.ones(int(Nt - 1)) + sigma_Y * dZt Yt = np.insert(np.exp(np.cumsum(yg)), 0, 1) Xt = np.ones(Nt) * nu * bet Deltabar = np.zeros(Nt) @@ -26,7 +26,9 @@ def BuildUpCohortsMAIN(dZt, Nt, dt, rho, nu, Vbar, mu_Y, sigma_Y, bet, That): tau = dt * np.ones(1) reduction = np.exp(-nu * dt) for i in range(1, Nt): - Part = IntVec * np.exp(-(rho + 0.5 * Delta_s_t * Delta_s_t) * dt + Delta_s_t * dZt[i - 1]) + Part = IntVec * np.exp( + -(rho + 0.5 * Delta_s_t * Delta_s_t) * dt + Delta_s_t * dZt[i - 1] + ) if i == 1: Xt[i] = Part Deltabar[i] = Part * Delta_s_t @@ -36,13 +38,14 @@ def BuildUpCohortsMAIN(dZt, Nt, dt, rho, nu, Vbar, mu_Y, sigma_Y, bet, That): IntVec = reduction * Part IntVec = np.append(IntVec, bet * (1 - reduction) * Xt[i]) f = IntVec / Xt[i] - dDelta_s_t = (PostVar(sigma_Y, Vbar, tau) / sigma_Y ** 2) * ( - -Delta_s_t * dt + np.ones(len(Delta_s_t)) * dZt[i - 1]) + dDelta_s_t = (PostVar(sigma_Y, Vbar, tau) / sigma_Y**2) * ( + -Delta_s_t * dt + np.ones(len(Delta_s_t)) * dZt[i - 1] + ) if i < Npre: Delta_s_t = Delta_s_t + dDelta_s_t Delta_s_t = np.append(Delta_s_t, 0) else: - DELbias = np.sum(dZt[int(i - Npre):i]) / That + DELbias = np.sum(dZt[int(i - Npre) : i]) / That Delta_s_t = Delta_s_t + dDelta_s_t Delta_s_t = np.append(Delta_s_t, DELbias) tau = tau + dt @@ -50,8 +53,24 @@ def BuildUpCohortsMAIN(dZt, Nt, dt, rho, nu, Vbar, mu_Y, sigma_Y, bet, That): return Deltabar, IntVec, Xt, Delta_s_t, Yt, Zt, f, tau -def SimCohortsMAIN(biasvec, dZt, Nt, tau, IntVec, Delta_s_t, dt, rho, nu, Vbar, mu_Y, sigma_Y, sigma_S, bet, That, - Npre): +def SimCohortsMAIN( + biasvec, + dZt, + Nt, + tau, + IntVec, + Delta_s_t, + dt, + rho, + nu, + Vbar, + mu_Y, + sigma_Y, + sigma_S, + bet, + That, + Npre, +): # Initializing variables Xt2 = np.ones(Nt) Deltabar2 = np.ones(Nt) @@ -66,7 +85,7 @@ def SimCohortsMAIN(biasvec, dZt, Nt, tau, IntVec, Delta_s_t, dt, rho, nu, Vbar, BIGF = np.zeros((Nt, Nt)) RevNt = np.flip(range(Nt)) - fhat = reduction ** RevNt * nu * dt + fhat = reduction**RevNt * nu * dt fhat = fhat / np.sum(fhat) # Expected returns @@ -82,18 +101,20 @@ def SimCohortsMAIN(biasvec, dZt, Nt, tau, IntVec, Delta_s_t, dt, rho, nu, Vbar, theta_t = np.zeros(Nt) fst = np.zeros(Nt) # individual connp.sumption-wealth ratio for i in range(Nt): - Part = IntVec * np.exp(-(0.5 * Delta_s_t ** 2) * dt + Delta_s_t * dZt[i]) + Part = IntVec * np.exp(-(0.5 * Delta_s_t**2) * dt + Delta_s_t * dZt[i]) if i < Nt - 1: if i == 0: part1[i] = np.sum(biasvec) / That else: - part1[i] = part1[i - 1] + (PostVar(sigma_Y, Vbar, (i * dt)) / sigma_Y ** 2) * ( - - part1[i - 1] * dt + dZt[i - 1]) + part1[i] = part1[i - 1] + ( + PostVar(sigma_Y, Vbar, (i * dt)) / sigma_Y**2 + ) * (-part1[i - 1] * dt + dZt[i - 1]) if i == 0: dR[i] = 0 else: - dR[i] = (mu_S[i - 1] - r_t[i - 1] + rho + mu_Y - sigma_Y ** 2 + nu * (1 - bet)) * dt + sigma_S * \ - dZt[i] + dR[i] = ( + mu_S[i - 1] - r_t[i - 1] + rho + mu_Y - sigma_Y**2 + nu * (1 - bet) + ) * dt + sigma_S * dZt[i] Xt2[i] = np.sum(Part) Deltabar2[i] = np.sum(Part * Delta_s_t) / Xt2[i] @@ -102,38 +123,65 @@ def SimCohortsMAIN(biasvec, dZt, Nt, tau, IntVec, Delta_s_t, dt, rho, nu, Vbar, BIGDELTA[i, :] = Delta_s_t mu_S[i] = sigma_S * sigma_Y - (sigma_S - sigma_Y) * Deltabar2[i] - mu_S_t[i] = sigma_S * sigma_Y + sigma_S * (-Deltabar2[i] + part1[i]) + sigma_Y * Deltabar2[i] + mu_S_t[i] = ( + sigma_S * sigma_Y + + sigma_S * (-Deltabar2[i] + part1[i]) + + sigma_Y * Deltabar2[i] + ) muhat_S_t[i] = mu_S[i] + sigma_S * np.sum(fhat * Delta_s_t) - r_t[i] = rho + mu_Y - sigma_Y ** 2 + nu * (1 - bet) + sigma_Y * Deltabar2[i] + r_t[i] = rho + mu_Y - sigma_Y**2 + nu * (1 - bet) + sigma_Y * Deltabar2[i] theta_t[i] = sigma_Y - Deltabar2[i] - Et[i] = np.sum(f * (Vbar / (1 + (Vbar / sigma_Y ** 2) * dt * RevNt)) * (1 / sigma_Y)) - Vt[i] = (np.sum(f * Delta_s_t ** 2) - Deltabar2[i] ** 2) * sigma_Y + Et[i] = np.sum( + f * (Vbar / (1 + (Vbar / sigma_Y**2) * dt * RevNt)) * (1 / sigma_Y) + ) + Vt[i] = (np.sum(f * Delta_s_t**2) - Deltabar2[i] ** 2) * sigma_Y fst[i] = f[Nt - i - 1] BIGF[i, :] = f - muC_s_t[i] = mu_Y + nu * (1 - bet) + (sigma_Y - Deltabar2[i]) * ( - part1[i] - Deltabar2[i]) + muC_s_t[i] = ( + mu_Y + nu * (1 - bet) + (sigma_Y - Deltabar2[i]) * (part1[i] - Deltabar2[i]) + ) sigmaC_s_t[i] = sigma_Y + part1[i] - Deltabar2[i] # Updating: - dDelta_s_t = (PostVar(sigma_Y, Vbar, tau) / sigma_Y ** 2) * ( - -Delta_s_t * dt + np.ones(len(Delta_s_t)) * dZt[i]) + dDelta_s_t = (PostVar(sigma_Y, Vbar, tau) / sigma_Y**2) * ( + -Delta_s_t * dt + np.ones(len(Delta_s_t)) * dZt[i] + ) if i < Npre: DELbias = (np.sum(biasvec[i:]) + np.sum(dZt[:i])) / That else: - DELbias = np.sum(dZt[i - Npre:i]) / That + DELbias = np.sum(dZt[i - Npre : i]) / That Delta_s_t = Delta_s_t[1:] + dDelta_s_t[1:] Delta_s_t = np.append(Delta_s_t, DELbias) IntVec = reduction * Part[1:] IntVec = np.append(IntVec, bet * (1 - reduction) * Xt2[i]) - Port = (part1 - Deltabar2) / sigma_S + (sigma_Y / sigma_S) * (1 - bet * np.flip(fhat) / fst) - - return Xt2, Deltabar2, part1, mu_S, mu_S_t, muhat_S_t, r_t, theta_t, Port, muC_s_t, sigmaC_s_t, BIGF, BIGDELTA, Et, Vt, dR + Port = (part1 - Deltabar2) / sigma_S + (sigma_Y / sigma_S) * ( + 1 - bet * np.flip(fhat) / fst + ) + + return ( + Xt2, + Deltabar2, + part1, + mu_S, + mu_S_t, + muhat_S_t, + r_t, + theta_t, + Port, + muC_s_t, + sigmaC_s_t, + BIGF, + BIGDELTA, + Et, + Vt, + dR, + ) ############################################################################################################# @@ -143,13 +191,15 @@ def SimCohortsMAIN(biasvec, dZt, Nt, tau, IntVec, Delta_s_t, dt, rho, nu, Vbar, nu = 0.02 # Death rate mu_Y = 0.02 # Growth rate of output sigma_Y = 0.033 # Standard deviation of output -sigma_S = sigma_Y # In equilibrium the stock price diffusion is the same as output diffusion +sigma_S = ( + sigma_Y # In equilibrium the stock price diffusion is the same as output diffusion +) w = 0.92 # Fraction of total output paid out as endowment # Some pre-calculations -D = rho ** 2 + 4 * (rho * nu + nu ** 2) * (1 - w) -bet = (rho + 2 * nu - D ** 0.5) / (2 * nu) -rlog = rho + mu_Y - sigma_Y ** 2 +D = rho**2 + 4 * (rho * nu + nu**2) * (1 - w) +bet = (rho + 2 * nu - D**0.5) / (2 * nu) +rlog = rho + mu_Y - sigma_Y**2 # Setting prior variance dt = 1 / 12 # time incremental @@ -161,7 +211,7 @@ def SimCohortsMAIN(biasvec, dZt, Nt, tau, IntVec, Delta_s_t, dt, rho, nu, Vbar, Nc = Nt T_hats = dt * np.arange(1, 13, 1) T_hat_dimension = len(T_hats) -#nus = [0.01, 0.02, 0.03] +# nus = [0.01, 0.02, 0.03] nus = [0.02] nu_dimension = len(nus) @@ -174,33 +224,60 @@ def SimCohortsMAIN(biasvec, dZt, Nt, tau, IntVec, Delta_s_t, dt, rho, nu, Vbar, f_parti_matrix = np.zeros((T_hat_dimension, nu_dimension, Mpaths)) for l in range(Mpaths): - dZ_build = dt ** 0.5 * np.random.randn(int(Nc - 1)) # dZt for the build function - dZ = dt ** 0.5 * np.random.randn(Nt) # dZt for the simulate function + dZ_build = dt**0.5 * np.random.randn(int(Nc - 1)) # dZt for the build function + dZ = dt**0.5 * np.random.randn(Nt) # dZt for the simulate function dZ_matrix[l, :] = dZ dZ_build_matrix[l, :] = dZ_build for k, That in enumerate(T_hats): - Vbar = (sigma_Y ** 2) / That # prior variance + Vbar = (sigma_Y**2) / That # prior variance Npre = int(That / dt) for m, nu in enumerate(nus): # this part is repetitive when there is only one value of nu - D = rho ** 2 + 4 * (rho * nu + nu ** 2) * (1 - w) - bet = (rho + 2 * nu - D ** 0.5) / (2 * nu) - rlog = rho + mu_Y - sigma_Y ** 2 + D = rho**2 + 4 * (rho * nu + nu**2) * (1 - w) + bet = (rho + 2 * nu - D**0.5) / (2 * nu) + rlog = rho + mu_Y - sigma_Y**2 # create age quartiles for analysis - ( - Deltabar, IntVec, Xt, Delta_s_t, Yt, Zt, f, tau - ) = BuildUpCohortsMAIN( + (Deltabar, IntVec, Xt, Delta_s_t, Yt, Zt, f, tau) = BuildUpCohortsMAIN( dZ_build, Nt, dt, rho, nu, Vbar, mu_Y, sigma_Y, bet, That ) dZforbias = np.diff(Zt) biasvec = dZforbias[-Npre:] - dZt = dt ** 0.5 * np.random.normal(0, 1, Nt) + dZt = dt**0.5 * np.random.normal(0, 1, Nt) Zt = np.cumsum(dZt) ( - Xt2, Deltabar2, Part1, mu_S, mu_S_t, muhat_S_t, r_t, theta_t, Port, muC_s_t, sigmaC_s_t, BIGf, BIGDELTA, - Et, Vt, dR + Xt2, + Deltabar2, + Part1, + mu_S, + mu_S_t, + muhat_S_t, + r_t, + theta_t, + Port, + muC_s_t, + sigmaC_s_t, + BIGf, + BIGDELTA, + Et, + Vt, + dR, ) = SimCohortsMAIN( - biasvec, dZt, Nt, tau, IntVec, Delta_s_t, dt, rho, nu, Vbar, mu_Y, sigma_Y, sigma_S, bet, That, Npre + biasvec, + dZt, + Nt, + tau, + IntVec, + Delta_s_t, + dt, + rho, + nu, + Vbar, + mu_Y, + sigma_Y, + sigma_S, + bet, + That, + Npre, ) r_matrix[k, m, l] = np.mean(r_t) theta_matrix[k, m, l] = np.mean(theta_t) @@ -211,7 +288,7 @@ def SimCohortsMAIN(biasvec, dZt, Nt, tau, IntVec, Delta_s_t, dt, rho, nu, Vbar, y1 = np.mean(r_matrix, axis=2) y2 = np.mean(theta_matrix, axis=2) -xlabels = ['interest rate', 'market price of risk'] +xlabels = ["interest rate", "market price of risk"] ys = [y1, y2] for i in range(len(ys)): @@ -220,8 +297,8 @@ def SimCohortsMAIN(biasvec, dZt, Nt, tau, IntVec, Delta_s_t, dt, rho, nu, Vbar, fig, ax = plt.subplots() # Create a figure containing a single axes. y = y_nu[:, j] ax.plot(x, y) - ax.set_xlabel('initial window') - ax.set_ylabel('mean ' + xlabels[i]) + ax.set_xlabel("initial window") + ax.set_ylabel("mean " + xlabels[i]) # plt.savefig('initial window and ' + xlabels[i] + '_' + mode + '_' + zoom_in + str(nu) + '.png', dpi=500, # format="png") # plt.savefig('initial window and ' + xlabels[i] + '_' + mode + '.png', dpi=500, format="png") diff --git a/src/cohort_builder.py b/src/cohort_builder.py index f637cbd..4c2fb6a 100644 --- a/src/cohort_builder.py +++ b/src/cohort_builder.py @@ -1,5 +1,11 @@ import numpy as np -from src.solver import bisection, solve_theta, find_the_rich, bisection_partial_constraint, solve_theta_partial_constraint +from src.solver import ( + bisection, + solve_theta, + find_the_rich, + bisection_partial_constraint, + solve_theta_partial_constraint, +) from tqdm import tqdm from typing import Tuple from src.stats import post_var, fadingmemo @@ -23,14 +29,7 @@ def build_cohorts( good_time_build: np.ndarray, mode_trade: str, mode_learn: str, -) -> Tuple[ - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, -]: +) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray,]: """builds up a sufficiently large set of cohorts in the economy, view each cohort as one agent with a constantly shrinking size Args: @@ -59,7 +58,7 @@ def build_cohorts( d_eta_st = np.zeros(1) # disagreement, eq(11) eta_bar = np.ones(1) eta_st_eta_ss = np.ones(1) - invest_tracker = np.ones(Ninit) if mode_trade == 'drop' else np.ones(Nc) + invest_tracker = np.ones(Ninit) if mode_trade == "drop" else np.ones(Nc) tau_info = np.ones(1) * dt for i in tqdm(range(1, Nc)): @@ -68,17 +67,22 @@ def build_cohorts( tau_short = tau[-i:] eta_st_eta_ss = eta_st_eta_ss * np.exp( - (-0.5 * d_eta_st ** 2) * dt - + d_eta_st * dZ_build[i - 1] + (-0.5 * d_eta_st**2) * dt + d_eta_st * dZ_build[i - 1] ) # equation (11) - eta_bar_parts = tax * np.exp(-tax * tau_short) * eta_bar * eta_st_eta_ss * dt # equation (18) - eta_bar_t = np.sum(eta_bar_parts) / ( 1 - tax * dt) # equation (18) # dividing by (1-tax*dt) keeps sum(f_st*dt) at 1 + eta_bar_parts = ( + tax * np.exp(-tax * tau_short) * eta_bar * eta_st_eta_ss * dt + ) # equation (18) + eta_bar_t = np.sum(eta_bar_parts) / ( + 1 - tax * dt + ) # equation (18) # dividing by (1-tax*dt) keeps sum(f_st*dt) at 1 # eta_bar_t = np.sum(eta_bar_parts) eta_st_eta_ss = np.append(eta_st_eta_ss, 1) eta_bar = np.append(eta_bar, eta_bar_t) - eta_bar = eta_bar / eta_bar_t # rescale, does not change the relative magnitude of each cohort + eta_bar = ( + eta_bar / eta_bar_t + ) # rescale, does not change the relative magnitude of each cohort # todo: eta_bar_t goes to 0 too quickly if (1) mode != 'comp', and (2) initial window very small # eta_bar_t is the denominator; it creates issues if too close to 0 # so I rescale eta_bar to keep it away from 0, without changing f_st @@ -87,59 +91,73 @@ def build_cohorts( f_st = np.append(f_st, tax) # update beliefs - dDelta_s_t = (post_var(sigma_Y, Vhat, tau_info) / sigma_Y**2 - ) * ( - -Delta_s_t * dt + dZ_build[i - 1] + dDelta_s_t = (post_var(sigma_Y, Vhat, tau_info) / sigma_Y**2) * ( + -Delta_s_t * dt + dZ_build[i - 1] ) # from eq(5) - if mode_learn == 'back_renew' and mode_trade == 'drop': + if mode_learn == "back_renew" and mode_trade == "drop": tau_info = np.append(tau_info, 0) + dt else: - tau_info = tau[-i - 1:] + tau_info = tau[-i - 1 :] if i < Npre: Delta_s_t += dDelta_s_t - Delta_s_t = np.append(Delta_s_t, 0) # newborns begin with 0 bias when there are not enough earlier observations + Delta_s_t = np.append( + Delta_s_t, 0 + ) # newborns begin with 0 bias when there are not enough earlier observations else: init_bias = np.average(dZ_build[int(i - Npre) : i]) / dt - Delta_s_t+= dDelta_s_t + Delta_s_t += dDelta_s_t Delta_s_t = np.append( Delta_s_t, init_bias ) # newborns begin with Npre earlier observations # find the market clearing theta, given beliefs and consumption shares - if i < Ninit or mode_trade == 'comp': # Ninit: initial rounds where the short-sale constraint is relaxed - d_eta_st = ( - Delta_s_t # relax the short-sale constraint in the beginning - ) + if ( + i < Ninit or mode_trade == "comp" + ): # Ninit: initial rounds where the short-sale constraint is relaxed + d_eta_st = Delta_s_t # relax the short-sale constraint in the beginning else: - invest_tracker = np.append(invest_tracker, 1) # all cohorts that are still in the market, new cohort by default can invest - if mode_trade == 'drop': + invest_tracker = np.append( + invest_tracker, 1 + ) # all cohorts that are still in the market, new cohort by default can invest + if mode_trade == "drop": if good_time_build[i - 1] == 1: - if mode_learn == 'back_collect': + if mode_learn == "back_collect": # agents who have left the market respond to the recent positive shocks # they collect all the information they missed during the drop period invest_tracker = np.ones(i + 1) # all can invest - if mode_learn == 'back_renew': - renew_bias = np.sum(dZ_build[int(i - Npre): i]) / (Npre * dt) - Delta_s_t = Delta_s_t * invest_tracker + renew_bias * (1 - invest_tracker) + if mode_learn == "back_renew": + renew_bias = np.sum(dZ_build[int(i - Npre) : i]) / (Npre * dt) + Delta_s_t = Delta_s_t * invest_tracker + renew_bias * ( + 1 - invest_tracker + ) tau_info = invest_tracker * tau_info + (1 - invest_tracker) * dt invest_tracker = np.ones(i + 1) possible_cons_share = f_st * dt * invest_tracker possible_delta_st = Delta_s_t * invest_tracker - lowest_bound = -np.max(possible_delta_st[np.nonzero(possible_delta_st)]) # absolute lower bound for theta among active investors + lowest_bound = -np.max( + possible_delta_st[np.nonzero(possible_delta_st)] + ) # absolute lower bound for theta among active investors theta_t = bisection( - solve_theta, lowest_bound, 50, possible_cons_share, possible_delta_st, sigma_Y + solve_theta, + lowest_bound, + 50, + possible_cons_share, + possible_delta_st, + sigma_Y, ) # solve for theta, 10 is a far away upper bound for theta a = Delta_s_t + theta_t - invest = (a >= 0) + invest = a >= 0 invest_tracker = invest * invest_tracker d_eta_st = a * invest_tracker - theta_t - if mode_trade == 'keep': - lowest_bound = -np.max(Delta_s_t[np.nonzero(Delta_s_t)]) # absolute lower bound for theta + if mode_trade == "keep": + lowest_bound = -np.max( + Delta_s_t[np.nonzero(Delta_s_t)] + ) # absolute lower bound for theta f_st_standard = f_st * dt theta_t = bisection( solve_theta, lowest_bound, 50, f_st_standard, Delta_s_t, sigma_Y @@ -159,27 +177,24 @@ def build_cohorts( ) - - - def build_cohorts_partial_constraint( - dZ_build: np.ndarray, - Nc: int, - dt: float, - tau: np.ndarray, - cohort_size: np.ndarray, - rho: float, - nu: float, - Vhat: float, - mu_Y: float, - sigma_Y: float, - tax: float, - Npre: int, - Ninit: int, - T_hat: float, - good_time_build: np.ndarray, - mode_trade: str, - mode_learn: str, + dZ_build: np.ndarray, + Nc: int, + dt: float, + tau: np.ndarray, + cohort_size: np.ndarray, + rho: float, + nu: float, + Vhat: float, + mu_Y: float, + sigma_Y: float, + tax: float, + Npre: int, + Ninit: int, + T_hat: float, + good_time_build: np.ndarray, + mode_trade: str, + mode_learn: str, ) -> Tuple[ np.ndarray, np.ndarray, @@ -224,25 +239,33 @@ def build_cohorts_partial_constraint( can_short_tracker = np.ones(Ninit) tau_info = np.ones(1) * dt - if mode_trade in ['rich_free'] and mode_learn in ['give_up', 'back_collect', 'back_renew']: + if mode_trade in ["rich_free"] and mode_learn in [ + "give_up", + "back_collect", + "back_renew", + ]: for i in tqdm(range(1, Nc)): # for i in tqdm(range(1, Ninit)): # new cohort born (age 0), get wealth transfer, observe, invest tau_short = tau[-i:] eta_st_eta_ss = eta_st_eta_ss * np.exp( - (-0.5 * d_eta_st ** 2) * dt - + d_eta_st * dZ_build[i - 1] + (-0.5 * d_eta_st**2) * dt + d_eta_st * dZ_build[i - 1] ) # equation (11) - eta_bar_parts = tax * np.exp(-tax * tau_short) * eta_bar * eta_st_eta_ss * dt # equation (18) + eta_bar_parts = ( + tax * np.exp(-tax * tau_short) * eta_bar * eta_st_eta_ss * dt + ) # equation (18) eta_bar_t = np.sum(eta_bar_parts) / ( - 1 - tax * dt) # equation (18) # dividing by (1-tax*dt) keeps sum(f_st*dt) at 1 + 1 - tax * dt + ) # equation (18) # dividing by (1-tax*dt) keeps sum(f_st*dt) at 1 # eta_bar_t = np.sum(eta_bar_parts) eta_st_eta_ss = np.append(eta_st_eta_ss, 1) eta_bar = np.append(eta_bar, eta_bar_t) - eta_bar = eta_bar / eta_bar_t # rescale, does not change the relative magnitude of each cohort + eta_bar = ( + eta_bar / eta_bar_t + ) # rescale, does not change the relative magnitude of each cohort # todo: eta_bar_t goes to 0 too quickly if (1) mode != 'comp', and (2) initial window very small # eta_bar_t is the denominator; it creates issues if too close to 0 # so I rescale eta_bar to keep it away from 0, without changing f_st @@ -252,22 +275,22 @@ def build_cohorts_partial_constraint( # update beliefs - dDelta_s_t = (post_var(sigma_Y, Vhat, tau_info) / sigma_Y ** 2 - ) * ( - -Delta_s_t * dt + dZ_build[i - 1] - ) # from eq(5) + dDelta_s_t = (post_var(sigma_Y, Vhat, tau_info) / sigma_Y**2) * ( + -Delta_s_t * dt + dZ_build[i - 1] + ) # from eq(5) - if mode_learn == 'back_renew': + if mode_learn == "back_renew": tau_info = np.append(tau_info, 0) + dt else: - tau_info = tau[-i - 1:] + tau_info = tau[-i - 1 :] if i < Npre: Delta_s_t += dDelta_s_t - Delta_s_t = np.append(Delta_s_t, - 0) # newborns begin with 0 bias when there are not enough earlier observations + Delta_s_t = np.append( + Delta_s_t, 0 + ) # newborns begin with 0 bias when there are not enough earlier observations else: - init_bias = np.average(dZ_build[int(i - Npre): i]) / dt + init_bias = np.average(dZ_build[int(i - Npre) : i]) / dt Delta_s_t += dDelta_s_t Delta_s_t = np.append( Delta_s_t, init_bias @@ -275,28 +298,30 @@ def build_cohorts_partial_constraint( # find the market clearing theta, given beliefs and consumption shares if i < Ninit: - d_eta_st = ( - Delta_s_t # relax the short-sale constraint in the beginning - ) + d_eta_st = Delta_s_t # relax the short-sale constraint in the beginning else: # add a new cohort in the trackers - invest_tracker = np.append(invest_tracker, - 1) # all cohorts that are still in the market, new cohort by default can invest - can_short_tracker = np.append(can_short_tracker, - 0) # some cohorts that are allowed to short, new cohort by default can't short - cohort_size_short = cohort_size[-i - 1:] + invest_tracker = np.append( + invest_tracker, 1 + ) # all cohorts that are still in the market, new cohort by default can invest + can_short_tracker = np.append( + can_short_tracker, 0 + ) # some cohorts that are allowed to short, new cohort by default can't short + cohort_size_short = cohort_size[-i - 1 :] if good_time_build[i - 1] == 1: - if mode_learn == 'back_collect': + if mode_learn == "back_collect": # agents who have left the market respond to the recent positive shocks # they collect all the information they missed during the drop period invest_tracker = np.ones(i + 1) # all can invest - if mode_learn == 'back_renew': - renew_bias = np.sum(dZ_build[int(i - Npre): i]) / (Npre * dt) - Delta_s_t = Delta_s_t * invest_tracker + renew_bias * (1 - invest_tracker) + if mode_learn == "back_renew": + renew_bias = np.sum(dZ_build[int(i - Npre) : i]) / (Npre * dt) + Delta_s_t = Delta_s_t * invest_tracker + renew_bias * ( + 1 - invest_tracker + ) tau_info = invest_tracker * tau_info + (1 - invest_tracker) * dt invest_tracker = np.ones(i + 1) @@ -307,26 +332,35 @@ def build_cohorts_partial_constraint( wealth_cutoff = find_the_rich( indiv_w_possible, cohort_size_possible, top=0.05 ) # find the cohorts that make the richest 1% pupolation in the current period that are still in the market - can_short = indiv_w_possible >= wealth_cutoff # these cohorts can short in this period + can_short = ( + indiv_w_possible >= wealth_cutoff + ) # these cohorts can short in this period can_short_tracker = ( - can_short_tracker + can_short >= 1 + can_short_tracker + can_short >= 1 ) # once rich, always can short theta_t = bisection_partial_constraint( - solve_theta_partial_constraint, -50, 50, can_short_tracker, Delta_s_t_possible, f_st_possible, - sigma_Y + solve_theta_partial_constraint, + -50, + 50, + can_short_tracker, + Delta_s_t_possible, + f_st_possible, + sigma_Y, ) want_to_short = (Delta_s_t + theta_t) < 0 - constrained = invest_tracker * want_to_short * ( - 1 - can_short_tracker + constrained = ( + invest_tracker * want_to_short * (1 - can_short_tracker) ) # in the market * want to short * can't short short_t = want_to_short * can_short_tracker long_t = (1 - want_to_short) * invest_tracker - invest_tracker = invest_tracker - constrained # constrained people drop, update invest_tracker + invest_tracker = ( + invest_tracker - constrained + ) # constrained people drop, update invest_tracker pi_st = (Delta_s_t + theta_t) / sigma_S * invest_tracker d_eta_st = Delta_s_t * invest_tracker - theta_t * (1 - invest_tracker) else: - print('mode not found') + print("mode not found") return ( Delta_s_t, @@ -339,8 +373,6 @@ def build_cohorts_partial_constraint( ) - - # # def build_cohorts_fading( # dZ_build: np.ndarray, @@ -471,4 +503,4 @@ def build_cohorts_partial_constraint( # intvec, # int_zt, # delta_ss, -# ) \ No newline at end of file +# ) diff --git a/src/cohort_simulator.py b/src/cohort_simulator.py index 84e6fa3..60740ab 100644 --- a/src/cohort_simulator.py +++ b/src/cohort_simulator.py @@ -1,39 +1,46 @@ import numpy as np from typing import Tuple from src.stats import post_var, fadingmemo -from src.solver import bisection, solve_theta, find_the_rich, bisection_partial_constraint, solve_theta_partial_constraint +from src.solver import ( + bisection, + solve_theta, + find_the_rich, + bisection_partial_constraint, + solve_theta_partial_constraint, +) from tqdm import tqdm from numba import jit + def simulate_cohorts( - Y: np.ndarray, - biasvec: np.ndarray, - dZ: np.ndarray, - Nt: int, - Nc: int, - tau: np.ndarray, - dt: float, - rho: float, - nu: float, - Vhat: float, - mu_Y: float, - sigma_Y: float, - sigma_S: float, - tax: float, - beta: float, - T_hat: float, - Npre: float, - Ninit: int, - mode_trade: str, - mode_learn: str, - cohort_size: np.ndarray, - Delta_s_t: np.ndarray, - eta_st_eta_ss: np.ndarray, - eta_bar: np.ndarray, - d_eta_st: np.ndarray, - invest_tracker: np.ndarray, - tau_info: np.ndarray, - good_time_simulate: np.ndarray, + Y: np.ndarray, + biasvec: np.ndarray, + dZ: np.ndarray, + Nt: int, + Nc: int, + tau: np.ndarray, + dt: float, + rho: float, + nu: float, + Vhat: float, + mu_Y: float, + sigma_Y: float, + sigma_S: float, + tax: float, + beta: float, + T_hat: float, + Npre: float, + Ninit: int, + mode_trade: str, + mode_learn: str, + cohort_size: np.ndarray, + Delta_s_t: np.ndarray, + eta_st_eta_ss: np.ndarray, + eta_bar: np.ndarray, + d_eta_st: np.ndarray, + invest_tracker: np.ndarray, + tau_info: np.ndarray, + good_time_simulate: np.ndarray, ) -> Tuple[ np.ndarray, np.ndarray, @@ -127,16 +134,19 @@ def simulate_cohorts( # new cohort born (age 0), get wealth transfer, observe, invest eta_st_eta_ss = eta_st_eta_ss * np.exp( - (-0.5 * d_eta_st ** 2) * dt - + d_eta_st * dZ_t + (-0.5 * d_eta_st**2) * dt + d_eta_st * dZ_t ) # equation (11) - eta_bar_parts = tax * np.exp(-tax * tau) * eta_bar * eta_st_eta_ss * dt # equation (18) + eta_bar_parts = ( + tax * np.exp(-tax * tau) * eta_bar * eta_st_eta_ss * dt + ) # equation (18) eta_bar_t = np.sum(eta_bar_parts) / (1 - tax * dt) # equation (18) eta_st_eta_ss = np.append(eta_st_eta_ss[1:], 1) eta_bar = np.append(eta_bar[1:], eta_bar_t) - eta_bar = eta_bar / eta_bar_t # rescale, does not change the relative magnitude of each cohort + eta_bar = ( + eta_bar / eta_bar_t + ) # rescale, does not change the relative magnitude of each cohort f_st = eta_bar_parts / eta_bar_t / dt f_st = np.append(f_st[1:], tax) @@ -146,59 +156,74 @@ def simulate_cohorts( w_cohort_st = Y[i] / beta * f_st w_st = w_cohort_st / cohort_size * dt else: - dR_t = mu_S_t * dt + sigma_S * dZ_t # realized stock return, mu_t^Sdt + sigma_t^Sdz_t + dR_t = ( + mu_S_t * dt + sigma_S * dZ_t + ) # realized stock return, mu_t^Sdt + sigma_t^Sdz_t w_t = Y[i] / beta - dw_st = ((r_t + nu - tax - beta) + pi_st * (mu_S_t - r_t)) * w_st * dt + w_st * pi_st * sigma_S * dZ_t # r_t, theta_t, pi_st from last loop, dZ_t just realized + dw_st = ( + (r_t + nu - tax - beta) + pi_st * (mu_S_t - r_t) + ) * w_st * dt + w_st * pi_st * sigma_S * dZ_t # r_t, theta_t, pi_st from last loop, dZ_t just realized w_st = w_st[1:] + dw_st[1:] w_st = np.append(w_st, w_t * tax / nu) w_cohort_st = w_st * cohort_size / dt # update beliefs - dDelta_s_t = ( - post_var(sigma_Y, Vhat, tau_info) / sigma_Y**2 - ) * ( + dDelta_s_t = (post_var(sigma_Y, Vhat, tau_info) / sigma_Y**2) * ( -Delta_s_t * dt + dZ_t ) # from eq(5) - if mode_learn == 'back_renew' and mode_trade == 'drop': + if mode_learn == "back_renew" and mode_trade == "drop": tau_info = np.append(tau_info[1:], 0) + dt - if i < Npre-1: - init_bias = (np.sum(biasvec[i+1:]) + np.sum(dZ[:i+1])) / T_hat + if i < Npre - 1: + init_bias = (np.sum(biasvec[i + 1 :]) + np.sum(dZ[: i + 1])) / T_hat else: - init_bias = np.sum(dZ[i+1 - Npre: i+1]) / T_hat + init_bias = np.sum(dZ[i + 1 - Npre : i + 1]) / T_hat Delta_s_t = Delta_s_t[1:] + dDelta_s_t[1:] Delta_s_t = np.append(Delta_s_t, init_bias) # find the market clearing theta, given beliefs and consumption shares of cohorts in the economy - invest_tracker = np.append(invest_tracker[1:], 1) # all cohorts that are still in the market, 1 by default + invest_tracker = np.append( + invest_tracker[1:], 1 + ) # all cohorts that are still in the market, 1 by default - if mode_trade == 'drop': + if mode_trade == "drop": if good_time_simulate[i] == 1: - if mode_learn == 'back_collect': + if mode_learn == "back_collect": # agents who have left the market respond to the recent positive shocks # they collect all the information they missed during the drop period invest_tracker = np.ones(Nc) # all can invest - if mode_learn == 'back_renew': + if mode_learn == "back_renew": if i < Npre - 1: - renew_bias = (np.sum(biasvec[i + 1:]) + np.sum(dZ[:i + 1])) / (Npre * dt) + renew_bias = ( + np.sum(biasvec[i + 1 :]) + np.sum(dZ[: i + 1]) + ) / (Npre * dt) else: - renew_bias = np.sum(dZ[i + 1 - Npre: i + 1]) / (Npre * dt) - Delta_s_t = Delta_s_t * invest_tracker + renew_bias * (1 - invest_tracker) + renew_bias = np.sum(dZ[i + 1 - Npre : i + 1]) / (Npre * dt) + Delta_s_t = Delta_s_t * invest_tracker + renew_bias * ( + 1 - invest_tracker + ) tau_info = invest_tracker * tau_info + (1 - invest_tracker) * dt invest_tracker = np.ones(Nc) # all can invest possible_cons_share = f_st * dt * invest_tracker possible_delta_st = Delta_s_t * invest_tracker - lowest_bound = -np.max(possible_delta_st) # absolute lower bound for theta among active investors + lowest_bound = -np.max( + possible_delta_st + ) # absolute lower bound for theta among active investors theta_t = bisection( - solve_theta, lowest_bound, 50, possible_cons_share, possible_delta_st, sigma_Y + solve_theta, + lowest_bound, + 50, + possible_cons_share, + possible_delta_st, + sigma_Y, ) # solve for theta a = Delta_s_t + theta_t - invest = (a >= 0) + invest = a >= 0 # want_to_short_t = invest_tracker * (1 - invest) invest_tracker = invest * invest_tracker d_eta_st = a * invest_tracker - theta_t @@ -210,15 +235,13 @@ def simulate_cohorts( age_t = np.sum(cohort_size * tau * invest_tracker) n_parti_t = np.sum(invest_tracker) / Nc - elif mode_trade == 'keep': + elif mode_trade == "keep": lowest_bound = -np.max(Delta_s_t) # absolute lower bound for theta f_st_standard = f_st * dt theta_t = bisection( solve_theta, lowest_bound, 50, f_st_standard, Delta_s_t, sigma_Y ) # solve for theta - d_eta_st = np.maximum( - -theta_t, Delta_s_t - ) # update max(Delta_s_t, -theta) + d_eta_st = np.maximum(-theta_t, Delta_s_t) # update max(Delta_s_t, -theta) invest = Delta_s_t >= -theta_t invest_fst = invest * f_st_standard popu_parti_t = np.sum(cohort_size * invest) @@ -228,7 +251,7 @@ def simulate_cohorts( age_t = np.sum(cohort_size * tau * invest) n_parti_t = np.sum(invest) / Nc - elif mode_trade == 'comp': + elif mode_trade == "comp": f_st_standard = f_st * dt Delta_bar_parti_t = np.sum(f_st_standard * Delta_s_t) theta_t = sigma_Y - Delta_bar_parti_t @@ -242,15 +265,10 @@ def simulate_cohorts( n_parti_t = np.sum(invest) / Nc else: - print('Warning! Mode trade not defined') + print("Warning! Mode trade not defined") exit() - r_t = ( - rho - + mu_Y - + nu - tax - - sigma_Y * theta_t - ) + r_t = rho + mu_Y + nu - tax - sigma_Y * theta_t mu_S_t = sigma_S * theta_t + r_t @@ -286,48 +304,36 @@ def simulate_cohorts( ) - - - - - - - - - - - - def simulate_cohorts_partial_constraint( - Y: np.ndarray, - biasvec: np.ndarray, - dZ: np.ndarray, - Nt: int, - Nc: int, - tau: np.ndarray, - dt: float, - rho: float, - nu: float, - Vhat: float, - mu_Y: float, - sigma_Y: float, - sigma_S: float, - tax: float, - beta: float, - T_hat: float, - Npre: float, - Ninit: int, - mode_trade: str, - mode_learn: str, - cohort_size: np.ndarray, - Delta_s_t: np.ndarray, - eta_st_eta_ss: np.ndarray, - eta_bar: np.ndarray, - d_eta_st: np.ndarray, - invest_tracker_t: np.ndarray, - can_short_tracker_t: np.ndarray, - tau_info: np.ndarray, - good_time_simulate: np.ndarray, + Y: np.ndarray, + biasvec: np.ndarray, + dZ: np.ndarray, + Nt: int, + Nc: int, + tau: np.ndarray, + dt: float, + rho: float, + nu: float, + Vhat: float, + mu_Y: float, + sigma_Y: float, + sigma_S: float, + tax: float, + beta: float, + T_hat: float, + Npre: float, + Ninit: int, + mode_trade: str, + mode_learn: str, + cohort_size: np.ndarray, + Delta_s_t: np.ndarray, + eta_st_eta_ss: np.ndarray, + eta_bar: np.ndarray, + d_eta_st: np.ndarray, + invest_tracker_t: np.ndarray, + can_short_tracker_t: np.ndarray, + tau_info: np.ndarray, + good_time_simulate: np.ndarray, ) -> Tuple[ np.ndarray, np.ndarray, @@ -447,16 +453,19 @@ def simulate_cohorts_partial_constraint( # new cohort born (age 0), get wealth transfer, observe, invest eta_st_eta_ss = eta_st_eta_ss * np.exp( - (-0.5 * d_eta_st ** 2) * dt - + d_eta_st * dZ_t + (-0.5 * d_eta_st**2) * dt + d_eta_st * dZ_t ) # equation (11) - eta_bar_parts = tax * np.exp(-tax * tau) * eta_bar * eta_st_eta_ss * dt # equation (18) + eta_bar_parts = ( + tax * np.exp(-tax * tau) * eta_bar * eta_st_eta_ss * dt + ) # equation (18) eta_bar_t = np.sum(eta_bar_parts) / (1 - tax * dt) # equation (18) eta_st_eta_ss = np.append(eta_st_eta_ss[1:], 1) eta_bar = np.append(eta_bar[1:], eta_bar_t) - eta_bar = eta_bar / eta_bar_t # rescale, does not change the relative magnitude of each cohort + eta_bar = ( + eta_bar / eta_bar_t + ) # rescale, does not change the relative magnitude of each cohort f_st = eta_bar_parts / eta_bar_t / dt f_st = np.append(f_st[1:], tax) @@ -479,46 +488,56 @@ def simulate_cohorts_partial_constraint( w_cohort_st = Y[i] / beta * f_st w_st = w_cohort_st / cohort_size * dt else: - dR_t = mu_S_t * dt + sigma_S * dZ_t # realized stock return, mu_t^Sdt + sigma_t^Sdz_t + dR_t = ( + mu_S_t * dt + sigma_S * dZ_t + ) # realized stock return, mu_t^Sdt + sigma_t^Sdz_t w_t = Y[i] / beta - dw_st = ((r_t + nu - tax - beta) + pi_st * (mu_S_t - r_t)) * w_st * dt + w_st * pi_st * sigma_S * dZ_t # r_t, theta_t, pi_st from last loop, dZ_t just realized + dw_st = ( + (r_t + nu - tax - beta) + pi_st * (mu_S_t - r_t) + ) * w_st * dt + w_st * pi_st * sigma_S * dZ_t # r_t, theta_t, pi_st from last loop, dZ_t just realized w_st = w_st[1:] + dw_st[1:] w_st = np.append(w_st, w_t * tax / nu) w_cohort_st = w_st * cohort_size / dt # update beliefs - dDelta_s_t = ( - post_var(sigma_Y, Vhat, tau_info) / sigma_Y**2 - ) * ( + dDelta_s_t = (post_var(sigma_Y, Vhat, tau_info) / sigma_Y**2) * ( -Delta_s_t * dt + dZ_t ) # from eq(5) - if mode == 'back_renew': + if mode == "back_renew": tau_info = np.append(tau_info[1:], 0) + dt if i < Npre - 1: - init_bias = (np.sum(biasvec[i+1:]) + np.sum(dZ[:i+1])) / T_hat + init_bias = (np.sum(biasvec[i + 1 :]) + np.sum(dZ[: i + 1])) / T_hat else: - init_bias = np.sum(dZ[i+1 - Npre: i+1]) / T_hat + init_bias = np.sum(dZ[i + 1 - Npre : i + 1]) / T_hat Delta_s_t = Delta_s_t[1:] + dDelta_s_t[1:] Delta_s_t = np.append(Delta_s_t, init_bias) - invest_tracker_t = np.append(invest_tracker_t[1:], 1) # all cohorts that are still in the market, 1 by default - can_short_tracker_t = np.append(can_short_tracker_t[1:], 0) # the cohorts that are allowed to short, 0 by default + invest_tracker_t = np.append( + invest_tracker_t[1:], 1 + ) # all cohorts that are still in the market, 1 by default + can_short_tracker_t = np.append( + can_short_tracker_t[1:], 0 + ) # the cohorts that are allowed to short, 0 by default if good_time_simulate[i] == 1: - if mode == 'back_collect': + if mode == "back_collect": # agents who have left the market respond to the recent positive shocks # they collect all the information they missed during the drop period invest_tracker_t = np.ones(Nc) # all can invest - if mode == 'back_renew': + if mode == "back_renew": if i < Npre - 1: - renew_bias = (np.sum(biasvec[i + 1:]) + np.sum(dZ[:i + 1])) / (Npre * dt) + renew_bias = (np.sum(biasvec[i + 1 :]) + np.sum(dZ[: i + 1])) / ( + Npre * dt + ) else: - renew_bias = np.sum(dZ[i + 1 - Npre: i + 1]) / (Npre * dt) - Delta_s_t = Delta_s_t * invest_tracker_t + renew_bias * (1 - invest_tracker_t) + renew_bias = np.sum(dZ[i + 1 - Npre : i + 1]) / (Npre * dt) + Delta_s_t = Delta_s_t * invest_tracker_t + renew_bias * ( + 1 - invest_tracker_t + ) tau_info = invest_tracker_t * tau_info + (1 - invest_tracker_t) * dt invest_tracker_t = np.ones(Nc) # all can invest @@ -530,30 +549,35 @@ def simulate_cohorts_partial_constraint( wealth_cutoff = find_the_rich( indiv_w_possible, cohort_size_possible, top=0.05 ) # find the cohorts that make the richest 1% pupolation in the current period that are still in the market - can_short = indiv_w_possible >= wealth_cutoff # these cohorts can short in this period + can_short = ( + indiv_w_possible >= wealth_cutoff + ) # these cohorts can short in this period can_short_tracker_t = ( can_short_tracker_t + can_short >= 1 ) # once rich, always can short theta_t = bisection_partial_constraint( - solve_theta_partial_constraint, -50, 50, can_short_tracker_t, Delta_s_t_possible, f_st_possible, sigma_Y + solve_theta_partial_constraint, + -50, + 50, + can_short_tracker_t, + Delta_s_t_possible, + f_st_possible, + sigma_Y, ) want_to_short = (Delta_s_t + theta_t) < 0 - constrained = invest_tracker_t * want_to_short * ( - 1 - can_short_tracker_t + constrained = ( + invest_tracker_t * want_to_short * (1 - can_short_tracker_t) ) # in the market * want to short * can't short short_t = want_to_short * can_short_tracker_t long_t = (1 - want_to_short) * invest_tracker_t - invest_tracker_t = invest_tracker_t - constrained # constrained people drop, update invest_tracker + invest_tracker_t = ( + invest_tracker_t - constrained + ) # constrained people drop, update invest_tracker pi_st = (Delta_s_t + theta_t) / sigma_S * invest_tracker_t d_eta_st = Delta_s_t * invest_tracker_t - theta_t * (1 - invest_tracker_t) - r_t = ( - rho - + mu_Y - + nu - tax - - sigma_Y * theta_t - ) + r_t = rho + mu_Y + nu - tax - sigma_Y * theta_t mu_S_t = sigma_S * theta_t + r_t @@ -570,7 +594,9 @@ def simulate_cohorts_partial_constraint( short[i, :] = short_t w_cohort[i, :] = w_cohort_st pi[i, :] = pi_st - Delta_bar_parti[i] = np.sum(Delta_s_t * invest_tracker_t * f_st) / np.sum(invest_tracker_t * f_st) + Delta_bar_parti[i] = np.sum(Delta_s_t * invest_tracker_t * f_st) / np.sum( + invest_tracker_t * f_st + ) Delta_bar_long[i] = np.sum(Delta_s_t * long_t * f_st) / np.sum(long_t * f_st) total_c_short = np.sum(short_t * f_st) if total_c_short == 0: @@ -620,10 +646,6 @@ def simulate_cohorts_partial_constraint( ) - - - - # # def simulate_cohorts_fading( # Y: np.ndarray, @@ -870,4 +892,4 @@ def simulate_cohorts_partial_constraint( # w_cohort, # age, # n_parti, -# ) \ No newline at end of file +# ) diff --git a/src/param.py b/src/param.py index ee3f891..46d5b6b 100644 --- a/src/param.py +++ b/src/param.py @@ -7,7 +7,7 @@ # nu = 0.03 mu_Y = 0.02 # Growth rate of output sigma_Y = 0.033 # Standard deviation of output -sigma_Y_sqr = sigma_Y ** 2 +sigma_Y_sqr = sigma_Y**2 sigma_S = ( sigma_Y # In equilibrium the stock price diffusion is the same as output diffusion ) @@ -23,7 +23,7 @@ dt = 1 / 12 # time incremental T_hat = 20 # Pre-trading period Npre = int(T_hat / dt) -Vhat = (sigma_Y ** 2) / T_hat # prior variance +Vhat = (sigma_Y**2) / T_hat # prior variance Ninit = int(20 / dt) T_cohort = 500 # time horizon to keep track of cohorts Nt = int(T_cohort / dt) # number of periods @@ -116,4 +116,4 @@ Et_matrix = np.zeros((Mpaths, Nt)) Vt_matrix = np.zeros((Mpaths, Nt)) -dR_matrix = np.zeros((Mpaths, Nt)) \ No newline at end of file +dR_matrix = np.zeros((Mpaths, Nt)) diff --git a/src/simulation.py b/src/simulation.py index 1d1090e..8a2fb75 100644 --- a/src/simulation.py +++ b/src/simulation.py @@ -6,7 +6,6 @@ from src.stats import shocks, good_times - def simulate( mode_trade: str, mode_learn: str, @@ -29,21 +28,21 @@ def simulate( tau: np.ndarray, cohort_size: np.ndarray, ) -> Tuple[ - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, ]: - ''' A program that combines cohort_builder and cohort_simulator, and finishes one whole simulation path + """A program that combines cohort_builder and cohort_simulator, and finishes one whole simulation path :param mode: scenario of the function, see param for scenario names :param Nc: number of cohorts :param dt: time increment @@ -56,7 +55,7 @@ def simulate( :param Npre: pre-trading periods :param T_hat: pre-trading years :return: - ''' + """ biasvec = dZ_build[-Npre:] # dZt used in the build_cohorts function @@ -65,12 +64,9 @@ def simulate( mu_Y, sigma_Y, dt, - ) + ) - ( - good_time_build, - good_time_simulate, - ) = good_times( + (good_time_build, good_time_simulate,) = good_times( dZ_build, dZ, dt, @@ -104,7 +100,7 @@ def simulate( good_time_build, mode_trade, mode_learn, - ) + ) ( r, @@ -149,7 +145,7 @@ def simulate( invest_tracker, tau_info_build, good_time_simulate, - ) + ) return ( r, @@ -190,34 +186,34 @@ def simulate_partial_constraint( tau: np.ndarray, cohort_size: np.ndarray, ) -> Tuple[ - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, - np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, ]: - ''' + """ :param mode: scenario of the function, see param for scenario names :param Nc: number of cohorts :param Nt: number of periods @@ -236,7 +232,7 @@ def simulate_partial_constraint( :param tau: age of each cohort :param cohort_size: size of each cohort :return: - ''' + """ # dZ_build = dt ** 0.5 * np.random.randn(int(Nc - 1)) # dZt for the build function biasvec = dZ_build[-Npre:] # dZt used in the build_cohorts function @@ -247,12 +243,9 @@ def simulate_partial_constraint( mu_Y, sigma_Y, dt, - ) + ) - ( - good_time_build, - good_time_simulate, - ) = good_times( + (good_time_build, good_time_simulate,) = good_times( dZ_build, dZ, dt, @@ -270,7 +263,25 @@ def simulate_partial_constraint( invest_tracker_build, can_short_tracker_build, tau_info_build, - ) = build_cohorts_partial_constraint(dZ_build, Nc, dt, tau, cohort_size, rho, nu, Vhat, mu_Y, sigma_Y, tax, Npre, Ninit, T_hat, good_time_build, mode_trade, mode_learn) + ) = build_cohorts_partial_constraint( + dZ_build, + Nc, + dt, + tau, + cohort_size, + rho, + nu, + Vhat, + mu_Y, + sigma_Y, + tax, + Npre, + Ninit, + T_hat, + good_time_build, + mode_trade, + mode_learn, + ) ( r, @@ -359,5 +370,3 @@ def simulate_partial_constraint( Delta_bar_long, Delta_bar_short, ) - - diff --git a/src/solver.py b/src/solver.py index 545c059..5930f7a 100644 --- a/src/solver.py +++ b/src/solver.py @@ -5,13 +5,13 @@ @jit(nopython=True) def bisection( - optimfun: Callable[[float, np.ndarray, np.ndarray, float], np.float64], - xlow: np.float64, - xhigh: np.float64, - arg1: np.ndarray, - arg2: np.ndarray, - arg3: float, - eps: float = 1e-9, + optimfun: Callable[[float, np.ndarray, np.ndarray, float], np.float64], + xlow: np.float64, + xhigh: np.float64, + arg1: np.ndarray, + arg2: np.ndarray, + arg3: float, + eps: float = 1e-9, ) -> np.float64: """Bisection method to solve x (theta) @@ -52,10 +52,10 @@ def bisection( @jit(nopython=True) def solve_theta( - theta_guess: np.float64, - consumption_share: np.ndarray, - Delta_s_t: np.ndarray, - sigma_Y: float, + theta_guess: np.float64, + consumption_share: np.ndarray, + Delta_s_t: np.ndarray, + sigma_Y: float, ) -> np.float64: """RHS - LHS of the eq(24), used to iteratively solve theta @@ -69,7 +69,7 @@ def solve_theta( np.float64: RHS - LHS """ invest = ( - Delta_s_t >= -theta_guess + Delta_s_t >= -theta_guess ) # eq(10) and eq(11), invest if theta_s_t >= -theta, constrained if otherwise invest_consumption_share = invest * consumption_share Delta_bar_parti = np.sum( @@ -82,21 +82,21 @@ def solve_theta( diff = 10000 else: diff = ( - sigma_Y - Delta_bar_parti - ) / total_invest_c_share - theta_guess # RHS - LHS, equals to 0 if find the right theta + sigma_Y - Delta_bar_parti + ) / total_invest_c_share - theta_guess # RHS - LHS, equals to 0 if find the right theta return diff -@jit(nopython = True) + +@jit(nopython=True) def find_the_rich( - indiv_w: np.ndarray, - cohort_size: np.ndarray, - top: float = 0.05) -> np.float64: - ''' + indiv_w: np.ndarray, cohort_size: np.ndarray, top: float = 0.05 +) -> np.float64: + """ :param indiv_w (np.ndarray): individual wealth of the agents, shape (Nc,) :param cohort_size (np.ndarray): shape (Nc,) :param top (float): can short criteria :return: a cutoff individual wealth level above which agents are then able to short - ''' + """ wealth_rank = indiv_w.argsort() indiv_w_sorted = indiv_w[wealth_rank[::-1]] cohort_size_sorted = cohort_size[wealth_rank[::-1]] @@ -106,15 +106,15 @@ def find_the_rich( return wealth_cutoff -@jit(nopython = True) +@jit(nopython=True) def solve_theta_partial_constraint( - theta_guess: float, - unconstrained: np.ndarray, - Delta_s_t: np.ndarray, - consumption_share: np.ndarray, - sigma_Y: float, + theta_guess: float, + unconstrained: np.ndarray, + Delta_s_t: np.ndarray, + consumption_share: np.ndarray, + sigma_Y: float, ) -> np.float64: - ''' + """ solve for theta in the conditionally constrained case, with the goal of market clearing in the stock market :param theta_guess (float): any guess of theta :param unconstrained (np.ndarray): the cohorts that can short @@ -122,26 +122,32 @@ def solve_theta_partial_constraint( :param consumption_share (np.ndarray): shape (Nc,) :param sigma_Y (float): volatility of aggregate output :return: the distance to converge - ''' + """ constrained = 1 - unconstrained - pi_constrained = np.maximum(Delta_s_t + theta_guess, 0) # investment if a cohort can't short - part_constrained = np.sum(pi_constrained * consumption_share * constrained) # for those cohorts that can't short - part_unconstrained = np.sum((Delta_s_t + theta_guess) * consumption_share * unconstrained) # for those cohorts that can short + pi_constrained = np.maximum( + Delta_s_t + theta_guess, 0 + ) # investment if a cohort can't short + part_constrained = np.sum( + pi_constrained * consumption_share * constrained + ) # for those cohorts that can't short + part_unconstrained = np.sum( + (Delta_s_t + theta_guess) * consumption_share * unconstrained + ) # for those cohorts that can short diff = (part_constrained + part_unconstrained) / sigma_Y - 1 return diff -@jit(nopython = True) +@jit(nopython=True) def solve_theta_partial_constraint_2( - theta_guess: float, - unconstrained: np.ndarray, - Delta_s_t: np.ndarray, - consumption_share: np.ndarray, - sigma_Y: float, + theta_guess: float, + unconstrained: np.ndarray, + Delta_s_t: np.ndarray, + consumption_share: np.ndarray, + sigma_Y: float, ) -> np.float64: - ''' + """ solve for theta in the conditionally constrained case, with the goal of market clearing in the stock market :param theta_guess (float): any guess of theta :param unconstrained (np.ndarray): the cohorts that can short @@ -149,9 +155,11 @@ def solve_theta_partial_constraint_2( :param consumption_share (np.ndarray): shape (Nc,) :param sigma_Y (float): volatility of aggregate output :return: the distance to converge - ''' + """ constrained = 1 - unconstrained - invest = (Delta_s_t >= -theta_guess) * constrained + unconstrained # eq(10) and eq(11), invest if theta_s_t >= -theta, constrained if otherwise + invest = ( + Delta_s_t >= -theta_guess + ) * constrained + unconstrained # eq(10) and eq(11), invest if theta_s_t >= -theta, constrained if otherwise invest_consumption_share = invest * consumption_share Delta_bar_parti = np.sum( Delta_s_t * invest_consumption_share @@ -163,23 +171,22 @@ def solve_theta_partial_constraint_2( diff = 10000 else: diff = ( - sigma_Y - Delta_bar_parti - ) / total_invest_c_share - theta_guess # RHS - LHS, equals to 0 if find the right theta + sigma_Y - Delta_bar_parti + ) / total_invest_c_share - theta_guess # RHS - LHS, equals to 0 if find the right theta return diff - @jit(nopython=True) def bisection_partial_constraint( - optimfun: Callable[[float, np.ndarray, np.ndarray, np.ndarray, float], np.float64], - xlow: np.float64, - xhigh: np.float64, - arg1: np.ndarray, - arg2: np.ndarray, - arg3: np.ndarray, - arg4: float, - eps: float = 1e-6, + optimfun: Callable[[float, np.ndarray, np.ndarray, np.ndarray, float], np.float64], + xlow: np.float64, + xhigh: np.float64, + arg1: np.ndarray, + arg2: np.ndarray, + arg3: np.ndarray, + arg4: float, + eps: float = 1e-6, ) -> np.float64: """Bisection method to solve x (theta) diff --git a/src/stats.py b/src/stats.py index c314762..20197e3 100644 --- a/src/stats.py +++ b/src/stats.py @@ -15,51 +15,42 @@ def post_var(sigma_Y: float, V_hat: float, tau: np.ndarray) -> np.ndarray: Returns: np.ndarray: shape (T, ) """ - sigma_Y_sq = sigma_Y ** 2 + sigma_Y_sq = sigma_Y**2 V = sigma_Y_sq * V_hat / (sigma_Y_sq + V_hat * tau) return V @jit(nopython=True) -def shocks( - dZt: np.ndarray, - mu_Y: float, - sigma_Y: float, - dt: float -) -> np.ndarray: +def shocks(dZt: np.ndarray, mu_Y: float, sigma_Y: float, dt: float) -> np.ndarray: """Calculate Zt and Yt from shocks dZt -Args: - dZt (np.ndarray): shocks for each period - mu_Y (float): mu_Y in eq(1) - sigma_Y (float): sigma_Y in eq(1) - dt (float): per unit of time + Args: + dZt (np.ndarray): shocks for each period + mu_Y (float): mu_Y in eq(1) + sigma_Y (float): sigma_Y in eq(1) + dt (float): per unit of time -Returns: - Zt (np.float64): cumulated shocks over the period - Yt (np.float64): output over the period -""" + Returns: + Zt (np.float64): cumulated shocks over the period + Yt (np.float64): output over the period""" # Zt = np.cumsum(dZt) # cumulated shocks, Nc * 1 # Zt = np.insert(Zt, 0, 0) yg = ( - mu_Y - 0.5 * sigma_Y ** 2 - ) * dt + sigma_Y * dZt # output in log, (Nc - 1) *1, eq(1) + mu_Y - 0.5 * sigma_Y**2 + ) * dt + sigma_Y * dZt # output in log, (Nc - 1) *1, eq(1) Yt = np.exp(np.cumsum(yg)) # Yt = np.insert(Yt, 0, 1) return Yt @jit(nopython=True) -def tau_calculator( - dt: float, - T_cohort: int -) -> np.ndarray: - """ Calculate tau -Args: - dt (float): per unit of time - T_cohort (int): number of periods -returns: - tau +def tau_calculator(dt: float, T_cohort: int) -> np.ndarray: + """Calculate tau + Args: + dt (float): per unit of time + T_cohort (int): number of periods + returns: + tau """ tau = np.arange(T_cohort, 0, -dt) return tau @@ -67,17 +58,14 @@ def tau_calculator( @jit(nopython=True) def good_times( - dZt_build: np.ndarray, - dZt: np.ndarray, - dt: float, - Nt: int, - Nc: int, - window: int, - z: float, -) -> Tuple[ - np.ndarray, - np.ndarray, -]: + dZt_build: np.ndarray, + dZt: np.ndarray, + dt: float, + Nt: int, + Nc: int, + window: int, + z: float, +) -> Tuple[np.ndarray, np.ndarray,]: """ returns the indicator for good times when agents previously dropped out from the stock market might return :param dZt_build: @@ -94,12 +82,12 @@ def good_times( if j < window: cummu_dZt_build[j] = 0 else: - cummu_dZt_build[j] = np.sum(dZt_build[j + 1 - window: j + 1]) + cummu_dZt_build[j] = np.sum(dZt_build[j + 1 - window : j + 1]) for i in range(Nt): if i < window: - cummu_dZt[i] = np.sum(dZt_build[i + 1 - window:]) + np.sum(dZt[: i + 1]) + cummu_dZt[i] = np.sum(dZt_build[i + 1 - window :]) + np.sum(dZt[: i + 1]) else: - cummu_dZt[i] = np.sum(dZt[i + 1 - window: i + 1]) + cummu_dZt[i] = np.sum(dZt[i + 1 - window : i + 1]) sigma_cummu = (dt * window) ** 0.5 good_time_build = cummu_dZt_build >= z * sigma_cummu good_time_simulate = cummu_dZt >= z * sigma_cummu @@ -107,9 +95,7 @@ def good_times( def fadingmemo(v, tau, sigma_Y, V_hat, int_zt, delta_ss): - v_st = np.log(1-v) / ( - (1-v) ** tau - 1 - ) - coef = v_st / (v_st * sigma_Y ** 2 + V_hat) - delta_st = coef * (sigma_Y ** 2 * delta_ss + V_hat * int_zt) - return delta_st \ No newline at end of file + v_st = np.log(1 - v) / ((1 - v) ** tau - 1) + coef = v_st / (v_st * sigma_Y**2 + V_hat) + delta_st = coef * (sigma_Y**2 * delta_ss + V_hat * int_zt) + return delta_st diff --git a/src/updates.py b/src/updates.py index 31e370d..e5ec75e 100644 --- a/src/updates.py +++ b/src/updates.py @@ -5,47 +5,46 @@ Npre = int(T_hat / dt) -Vhat = (sigma_Y ** 2) / T_hat # prior variance +Vhat = (sigma_Y**2) / T_hat # prior variance DELbias_mat = np.zeros(Nc) Delta_s_t_mat = np.zeros((Nc, Nt)) -T_hatS = dt * np.arange(1,13,1) +T_hatS = dt * np.arange(1, 13, 1) + def updates(tau, sigma_Y, dZ_build, dZ, T_hat): Delta_s_t = np.zeros(1) Npre = int(T_hat / dt) - Vhat = (sigma_Y ** 2) / T_hat # prior variance + Vhat = (sigma_Y**2) / T_hat # prior variance DELbias_mat = np.zeros(Nt) Delta_s_t_mat = np.zeros((Nc, Nt)) biasvec = dZ_build[-Npre:] for i in range(1, Nc): tau_short = tau[-i] - dDelta_s_t = ( - post_var(sigma_Y, Vhat, tau_short) / sigma_Y ** 2 - ) * ( - -Delta_s_t * dt + dZ_build[i - 1] - ) # from eq(5) + dDelta_s_t = (post_var(sigma_Y, Vhat, tau_short) / sigma_Y**2) * ( + -Delta_s_t * dt + dZ_build[i - 1] + ) # from eq(5) if i < Npre: Delta_s_t += dDelta_s_t - Delta_s_t = np.append(Delta_s_t, 0) # newborns begin with 0 bias when there are not enough earlier observations + Delta_s_t = np.append( + Delta_s_t, 0 + ) # newborns begin with 0 bias when there are not enough earlier observations else: - DELbias = np.average(dZ_build[int(i - Npre): i]) / dt + DELbias = np.average(dZ_build[int(i - Npre) : i]) / dt Delta_s_t += dDelta_s_t Delta_s_t = np.append( Delta_s_t, DELbias ) # newborns begin with Npre earlier observations for j in range(Nt): - dDelta_s_t = ( - post_var(sigma_Y, Vhat, tau) / sigma_Y**2 - ) * ( + dDelta_s_t = (post_var(sigma_Y, Vhat, tau) / sigma_Y**2) * ( -Delta_s_t * dt + dZ[j] ) # from eq(5) if j < (Npre - 1): - init_bias = (np.sum(biasvec[j+1:]) + np.sum(dZ[:j+1])) / T_hat + init_bias = (np.sum(biasvec[j + 1 :]) + np.sum(dZ[: j + 1])) / T_hat else: - init_bias = np.average(dZ[j+1 - Npre: j+1]) / dt + init_bias = np.average(dZ[j + 1 - Npre : j + 1]) / dt DELbias_mat[j] = init_bias Delta_s_t += dDelta_s_t @@ -54,6 +53,7 @@ def updates(tau, sigma_Y, dZ_build, dZ, T_hat): return np.mean(Delta_s_t_mat[:, -100:]) + m = 200 n = len(T_hatS) @@ -61,8 +61,8 @@ def updates(tau, sigma_Y, dZ_build, dZ, T_hat): for a in range(m): - dZ_build = dt ** 0.5 * np.random.randn(int(Nc - 1)) - dZ = dt ** 0.5 * np.random.randn(Nt) + dZ_build = dt**0.5 * np.random.randn(int(Nc - 1)) + dZ = dt**0.5 * np.random.randn(Nt) for b, t_hat in enumerate(T_hatS): # print(t_hat) A[a, b] = updates(tau, sigma_Y, dZ_build, dZ, t_hat) diff --git a/tax rate experiment.py b/tax rate experiment.py index aaa1bb3..c0bbf94 100644 --- a/tax rate experiment.py +++ b/tax rate experiment.py @@ -8,7 +8,7 @@ # modes = ['drop'] -modes = ['rich_free'] +modes = ["rich_free"] tax_rates = np.arange(0.005, 0.021, 0.002) tax_rate_dimension = len(tax_rates) @@ -62,8 +62,10 @@ for l in range(Mpaths): s = time.time() # same shocks for the different T_hats - dZ_build = dt ** 0.5 * np.random.randn(int(Nc - 1)) # dZt for the build function - dZ = dt ** 0.5 * np.random.randn(Nt) # dZt for the simulate function + dZ_build = dt**0.5 * np.random.randn( + int(Nc - 1) + ) # dZt for the build function + dZ = dt**0.5 * np.random.randn(Nt) # dZt for the simulate function dZ_matrix[l, :] = dZ Z = np.cumsum(dZ) Z_matrix[l, :] = Z @@ -71,7 +73,7 @@ for k, tax_rate in enumerate(tax_rates): omega = rho + nu - tax_rate - if mode == 'drop': + if mode == "drop": ( mu_S, mu_S_s, @@ -89,8 +91,25 @@ w_cohort, age_parti, n_parti, - ) = simulate(mode, Nc, Nt, dt, rho, nu, Vhat, mu_Y, sigma_Y, tax_rate, omega, Npre, T_hat, dZ_build, dZ, tau, - cohort_size) + ) = simulate( + mode, + Nc, + Nt, + dt, + rho, + nu, + Vhat, + mu_Y, + sigma_Y, + tax_rate, + omega, + Npre, + T_hat, + dZ_build, + dZ, + tau, + cohort_size, + ) dR_matrix[k, l] = np.average(dR[1200:]) r_matrix[k, l] = np.average(r[1200:]) @@ -98,9 +117,15 @@ popu_parti_matrix[k, l] = np.average(popu_parti[1200:]) invest = pi > 0 parti_rate = invest * cohort_size - popu_age1_matrix[k, l] = np.average(np.sum(parti_rate[1200:, tau_cutoff1:], axis=1)) - popu_age2_matrix[k, l] = np.average(np.sum(parti_rate[1200:, tau_cutoff2:], axis=1)) - popu_age3_matrix[k, l] = np.average(np.sum(parti_rate[1200:, tau_cutoff3:], axis=1)) + popu_age1_matrix[k, l] = np.average( + np.sum(parti_rate[1200:, tau_cutoff1:], axis=1) + ) + popu_age2_matrix[k, l] = np.average( + np.sum(parti_rate[1200:, tau_cutoff2:], axis=1) + ) + popu_age3_matrix[k, l] = np.average( + np.sum(parti_rate[1200:, tau_cutoff3:], axis=1) + ) age_parti_matrix[k, l] = np.average(age_parti[1200:]) n_parti_matrix[k, l] = np.average(n_parti[1200:]) @@ -135,8 +160,25 @@ Delta_bar_parti, Delta_bar_long, Delta_bar_short, - ) = simulate_partial_constraint(mode, Nc, Nt, dt, rho, nu, Vhat, mu_Y, sigma_Y, tax_rate, omega, Npre, - T_hat, dZ_build, dZ, tau, cohort_size) + ) = simulate_partial_constraint( + mode, + Nc, + Nt, + dt, + rho, + nu, + Vhat, + mu_Y, + sigma_Y, + tax_rate, + omega, + Npre, + T_hat, + dZ_build, + dZ, + tau, + cohort_size, + ) dR_matrix[k, l] = np.average(dR[1200:]) r_matrix[k, l] = np.average(r[1200:]) @@ -144,15 +186,20 @@ popu_parti_matrix[k, l] = np.average(popu_parti[1200:]) invest = pi > 0 parti_rate = invest * cohort_size - popu_age1_matrix[k, l] = np.average(np.sum(parti_rate[1200:, tau_cutoff1:], axis=1)) - popu_age2_matrix[k, l] = np.average(np.sum(parti_rate[1200:, tau_cutoff2:], axis=1)) - popu_age3_matrix[k, l] = np.average(np.sum(parti_rate[1200:, tau_cutoff3:], axis=1)) + popu_age1_matrix[k, l] = np.average( + np.sum(parti_rate[1200:, tau_cutoff1:], axis=1) + ) + popu_age2_matrix[k, l] = np.average( + np.sum(parti_rate[1200:, tau_cutoff2:], axis=1) + ) + popu_age3_matrix[k, l] = np.average( + np.sum(parti_rate[1200:, tau_cutoff3:], axis=1) + ) age_parti_matrix[k, l] = np.average(age_parti[1200:]) n_parti_matrix[k, l] = np.average(n_parti[1200:]) # covariance: - # graphs: x = tax_rates # y0 = (np.ones(tax_rate_dimension) * sigma_Y ** 2) / x @@ -166,33 +213,47 @@ y5 = np.average(n_parti_matrix, axis=1) y6 = -y2 * sigma_Y + mu_Y - xlabels = ['interest rate', 'market price of risk', 'participation rate', 'age of participants', - 'number of cohorts', 'cutoff belief'] + xlabels = [ + "interest rate", + "market price of risk", + "participation rate", + "age of participants", + "number of cohorts", + "cutoff belief", + ] ys = [y1, y2, y3, y4, y5, y6] for i in range(len(ys)): fig, ax = plt.subplots() # Create a figure containing a single axes. y = ys[i] if i == 2: - ax.fill_between(x, y31, color='steelblue', linewidth=0.4, label='20 < Age <= 35, youngest quartile') - ax.fill_between(x, y32, y31, color='darkseagreen', linewidth=0.4, label='35 < Age <= 55') - ax.fill_between(x, y33, y32, color='moccasin', linewidth=0.4, label='55 < Age <= 89') - ax.fill_between(x, y, y33, color='pink', linewidth=0.4, label='Age > 89, oldest quartile') + ax.fill_between( + x, + y31, + color="steelblue", + linewidth=0.4, + label="20 < Age <= 35, youngest quartile", + ) + ax.fill_between( + x, y32, y31, color="darkseagreen", linewidth=0.4, label="35 < Age <= 55" + ) + ax.fill_between( + x, y33, y32, color="moccasin", linewidth=0.4, label="55 < Age <= 89" + ) + ax.fill_between( + x, + y, + y33, + color="pink", + linewidth=0.4, + label="Age > 89, oldest quartile", + ) else: ax.plot(x, y) - ax.set_xlabel('tax rate') + ax.set_xlabel("tax rate") if i == 5: ax.set_ylabel(xlabels[i]) else: - ax.set_ylabel('mean ' + xlabels[i]) - plt.savefig('tax rate and ' + xlabels[i] + '_' + mode + '.png', dpi=500) + ax.set_ylabel("mean " + xlabels[i]) + plt.savefig("tax rate and " + xlabels[i] + "_" + mode + ".png", dpi=500) # plt.savefig('initial window and ' + xlabels[i] + '_' + mode + '.png', dpi=500) - - - - - - - - -