<jupyter_start><jupyter_text>Monte Carlo Simulation of the Nonparametric and Semiparametric Estimators<jupyter_code># load general packages
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
import numba as nb

# Import user package
import functions as funcs

## Uncomment for nicer plots ##
plt.rc('font', family='serif')
plt.rc('text', usetex=True)
plt.rc('axes', axisbelow=True)
# Set seed once
np.seed = 200
# Set Monte Carlo parameters
reps = 10_000 # Repetitions
n = 50_000 # Size of each sample<jupyter_output><empty_output><jupyter_text>For $\alpha = 1.5$<jupyter_code># Set distribution parameters
alpha_1 = 1.5
mu = 0.0

# Initialize
G_vec_1 = np.zeros(shape=reps)
G_nonparam_vec_1 = np.zeros(shape=reps)

# Loop over replications for now
for rep in range(reps):
    
    # Draw values
    rand = np.random.uniform(0,1, size=n)
    x = mu - 1 + (1 - rand)**(-1/alpha_1)
    
    # Compute p
    P = 0.01
    u = np.percentile(x,100*(1-P))

    # Estmate alpha
    N_1 = np.size(x[x>=u])    
    alpha_hat = N_1 / (np.sum(np.log(x[x>u])) - N_1*np.log(u))
    
    # Compute S
    x_low = x[x<u] # Is this needed?
    E_low = np.mean(x[x<u])
    S = alpha_hat/(alpha_hat-1)*u*P / (alpha_hat/(alpha_hat-1)*u*P + E_low * (1-P))
    S_alvaredo = np.mean(x[x>=u]) / (np.mean(x[x>=u]) + E_low) # Numerical top income share

    # Compute the Gini coefficients with both methods
    G_high = 1/(2*alpha_hat-1)
    G_low = funcs.gini(x[x<u])
    
    G = G_high*P*S + G_low*(1-P)*(1-S) + S - P
    G_nonparam = funcs.gini(x)
    
    # Save values
    G_vec_1[rep] = G
    G_nonparam_vec_1[rep] = G_nonparam<jupyter_output><empty_output><jupyter_text>Plot the results<jupyter_code># Compute the true gini
eta_1 = mu-1 + alpha_1/(alpha_1-1)
G_true_1 = 1 - 2/eta_1 * ( (mu-1)/2 + alpha_1/(alpha_1-1)*(1 - alpha_1/(2*alpha_1-1)))

# Plot
fig = plt.figure(figsize=(6,4))
ax = fig.add_subplot(1,1,1)

ax.grid(b=True, which = 'major', linestyle='-', linewidth=0.5, color='0.9', zorder=0)
ax.tick_params(axis='both', bottom=True, top=True, left=True, right=True, direction='in', which='both')

ax.hist(G_vec_1, density = True, bins = np.linspace(.7,.85,200), color = 'red', alpha = 0.7, zorder = 2, label='Semiparametric', histtype='stepfilled')
ax.hist(G_nonparam_vec_1, density = True, bins = np.linspace(.7,.85,200), color = 'blue', alpha = 0.7, zorder = 2, label='Nonparametric', histtype='stepfilled')

ax.vlines(G_true_1, ymin=0,ymax=100, linestyle = '-', linewidth = 1.0, label='True Value', zorder=3)
ax.set_ylim(0,60)
ax.set_xlim(0.7,.85)

ax.legend(frameon = True, edgecolor = 'k', facecolor = 'white', framealpha=1, fancybox=False)

plt.savefig('../../tex/figs/estimator_15.pdf')<jupyter_output><empty_output><jupyter_text>Compute kernel density<jupyter_code># First compute the mode using kernel density
kd_nonparam_1 = stats.gaussian_kde(G_nonparam_vec_1)
kd_semiparam_1 = stats.gaussian_kde(G_vec_1)

# Define point on which to evaluate the estimated pdf's
points_1 = np.linspace(.71,.85,200)
kd_11 = kd_nonparam_1.pdf(points)
kd_12 = kd_semiparam_1.pdf(points)
# Plot results
# Plot
fig = plt.figure(figsize=(6,4))
ax = fig.add_subplot(1,1,1)

ax.grid(b=True, which = 'major', linestyle='-', linewidth=0.5, color='0.9', zorder=0)
ax.tick_params(axis='both', bottom=True, top=True, left=True, right=True, direction='in', which='both')

ax.plot(points_1, kd_11, color = 'blue', linestyle ='-', zorder = 2, linewidth=1.5, label = 'Nonparametric')
ax.plot(points_1, kd_12, color = 'red', linestyle = '-', zorder = 2, linewidth = 1.5, label = 'Semiparametric')

ax.vlines(G_true_1, ymin=0,ymax=100, linestyle = '-', linewidth = 1.0, label='True Value', zorder=3)
ax.set_ylim(0,60)
ax.set_xlim(0.7,.85)

ax.legend(frameon = True, edgecolor = 'k', facecolor = 'white', framealpha=1, fancybox=False)

plt.savefig('../../tex/figs/kd_1.pdf')
# Compute mode
mode_nonparam_1 = points_1[np.argmax(kd_11)]
mode_semiparam_1 = points_1[np.argmax(kd_12)]<jupyter_output><empty_output><jupyter_text>Compute summary statistics<jupyter_code>print('Nonparam bias = ' + str(np.mean(G_nonparam_vec_1) - G_true_1))
print('Semiparam bias = ' + str(np.mean(G_vec_1) - G_true_1))

print('Nonparam bias (mode) = ' + str(mode_nonparam_1 - G_true_1))
print('Semiparam bias (mode) = ' + str(mode_semiparam_1 - G_true_1))

print('Nonparam variance = ' + str(np.std(G_nonparam_vec_1)))
print('Semiparam variance = ' + str(np.std(G_vec_1)))

print('Nonparam skewness = ' + str(stats.skew(G_nonparam_vec_1)))
print('Semiparam skewness = ' + str(stats.skew(G_vec_1)))<jupyter_output>Nonparam bias = -0.0018597624801071078
Semiparam bias = 0.002628326715806817
Nonparam bias (mode) = -0.0069346733668341765
Semiparam bias (mode) = 0.00010050251256277232
Nonparam variance = 0.01533912639731005
Semiparam variance = 0.008577228784523399
Nonparam skewness = 3.920432265624547
Semiparam skewness = 0.468590213803735<jupyter_text>**See if I can compute the mode somehow** For $\alpha = 1.2$<jupyter_code># Set distribution parameters
alpha_2 = 1.2
mu = 0.0

# Initialize
G_vec_2 = np.zeros(shape=reps)
G_nonparam_vec_2 = np.zeros(shape=reps)

# Loop over replications for now
for rep in range(reps):
    
    # Draw values
    rand = np.random.uniform(0,1, size=n)
    x = mu - 1 + (1 - rand)**(-1/alpha_2)
    
    # Compute p
    P = 0.01
    u = np.percentile(x,100*(1-P))

    # Estmate alpha
    N_1 = np.size(x[x>=u])    
    alpha_hat = N_1 / (np.sum(np.log(x[x>u])) - N_1*np.log(u))
    
    # Compute S
    x_low = x[x<u] # Is this needed?
    E_low = np.mean(x[x<u])
    S = alpha_hat/(alpha_hat-1)*u*P / (alpha_hat/(alpha_hat-1)*u*P + E_low * (1-P))
    S_alvaredo = np.mean(x[x>=u]) / (np.mean(x[x>=u]) + E_low) # Numerical top income share

    # Compute the Gini coefficients with both methods
    G_high = 1/(2*alpha_hat-1)
    G_low = funcs.gini(x[x<u])
    
    G = G_high*P*S + G_low*(1-P)*(1-S) + S - P
    G_nonparam = funcs.gini(x)
    
    # Save values
    G_vec_2[rep] = G
    G_nonparam_vec_2[rep] = G_nonparam<jupyter_output><empty_output><jupyter_text>Plot the results<jupyter_code># Compute the true gini
eta_2 = mu-1 + alpha_2/(alpha_2-1)
G_true_2 = 1 - 2/eta_2 * ( (mu-1)/2 + alpha_2/(alpha_2-1)*(1 - alpha_2/(2*alpha_2-1)))

# Plot
fig = plt.figure(figsize=(6,4))
ax = fig.add_subplot(1,1,1)

ax.grid(b=True, which = 'major', linestyle='-', linewidth=0.5, color='0.9', zorder=0)
ax.tick_params(axis='both', bottom=True, top=True, left=True, right=True, direction='in', which='both')

ax.hist(G_vec_2, density = True, bins = np.linspace(.75,1,200), color = 'red', alpha = 0.7, zorder = 2, label='Semiparametric', histtype='stepfilled')
ax.hist(G_nonparam_vec_2, density = True, bins =  np.linspace(.75,1,200), color = 'blue', alpha = 0.7, zorder = 2, label='Nonparametric', histtype='stepfilled')

ax.vlines(G_true_2, ymin=0,ymax=100, linestyle = '-', linewidth = 1.0, label='True Value', zorder=3)
ax.set_ylim(0,30)
ax.set_xlim(0.75,1)

ax.legend(frameon = True, edgecolor = 'k', facecolor = 'white', framealpha=1, fancybox=False)

plt.savefig('../../tex/figs/estimator_12.pdf')
# First compute the mode using kernel density
kd_nonparam_2 = stats.gaussian_kde(G_nonparam_vec_2)
kd_semiparam_2 = stats.gaussian_kde(G_vec_2)

# Define point on which to evaluate the estimated pdf's
points_2 = np.linspace(.75,1,200)
kd_21 = kd_nonparam_2.pdf(points_2)
kd_22 = kd_semiparam_2.pdf(points_2)
# Plot results
# Plot
fig = plt.figure(figsize=(6,4))
ax = fig.add_subplot(1,1,1)

ax.grid(b=True, which = 'major', linestyle='-', linewidth=0.5, color='0.9', zorder=0)
ax.tick_params(axis='both', bottom=True, top=True, left=True, right=True, direction='in', which='both')

ax.plot(points_2, kd_21, color = 'blue', linestyle ='-', zorder = 2, linewidth=1.5, label = 'Nonparametric')
ax.plot(points_2, kd_22, color = 'red', linestyle = '-', zorder = 2, linewidth = 1.5, label = 'Semiparametric')

ax.vlines(G_true_2, ymin=0,ymax=100, linestyle = '-', linewidth = 1.0, label='True Value', zorder=3)
ax.set_ylim(0,30)
ax.set_xlim(0.75,1)

ax.legend(frameon = True, edgecolor = 'k', facecolor = 'white', framealpha=1, fancybox=False)

plt.savefig('../../tex/figs/kd_2.pdf')
# Compute mode
mode_nonparam_2 = points_2[np.argmax(kd_21)]
mode_semiparam_2 = points_2[np.argmax(kd_22)]
print('Nonparam bias = ' + str(np.mean(G_nonparam_vec_2) - G_true_2))
print('Semiparam bias = ' + str(np.mean(G_vec_2) - G_true_2))

print('Nonparam bias (mode) = ' + str(mode_nonparam_2 - G_true_2))
print('Semiparam bias (mode) = ' + str(mode_semiparam_2 - G_true_2))

print('Nonparam variance = ' + str(np.std(G_nonparam_vec_2)))
print('Semiparam variance = ' + str(np.std(G_vec_2)))

print('Nonparam skewness = ' + str(stats.skew(G_nonparam_vec_2)))
print('Semiparam skewness = ' + str(stats.skew(G_vec_2)))<jupyter_output>Nonparam bias = -0.017410439265306055
Semiparam bias = 0.005315192861668483
Nonparam bias (mode) = -0.029253409906676464
Semiparam bias (mode) = 0.0008973438621677765
Nonparam variance = 0.026030803086283336
Semiparam variance = 0.020248307353577407
Nonparam skewness = 2.0481043843392492
Semiparam skewness = 0.7135023025704469
{'fraction_non_alphanumeric': 0.16575492341356673, 'fraction_numerical': 0.07330415754923414, 'mean_word_length': 4.837164750957855, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '352436', 'license': 'MIT', 'n_tokens_mistral': 4163, 'n_tokens_neox': 3842, 'n_words': 714, 'path': '01 - EstimatorMC.ipynb', 'repo_name': 'madsankern/TermPaper_ineq', 'size': '9140'}