<jupyter_start><jupyter_text>Importance Sampling and Monte Carlo for the Heterogenous Agent Model<jupyter_code># Load general packages
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import interpolate

# Import user generated functions
import functions as funcs

## Uncomment for nicer plots ##
plt.rc('font', family='serif')
plt.rc('text', usetex=True)
plt.rc('axes', axisbelow=True)<jupyter_output>The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload<jupyter_text>Load model data and plot Load the data and save the wealth distribution. The distributions is generated from an external matlab file 'model.m'.<jupyter_code># Load the matlab generated wealth distribution
df_g = pd.read_csv(r'ga_var.csv', header=None)
df_a = pd.read_csv(r'a_var.csv', header=None)

# Save as a numpy array
g_bottom = df_g.to_numpy()
w_bottom = df_a.to_numpy()
g_bottom = np.array(g_bottom).ravel() # Remove outer brackets
w_bottom = np.array(w_bottom).ravel()
# Define variables
alpha = 1.5
mu = -0.3
np.random.seed(214)<jupyter_output><empty_output><jupyter_text>Manually add a right tail to the wealth distribution<jupyter_code># Normalization constant for upper tail
C = g_bottom[-1] / 2000**(-alpha-1)

# Grid for tail x values
tail_vals = np.linspace(w_bottom[-1],1_000_000, 100000)

# Compute top
g_top = C*tail_vals**(-alpha-1) # alpha value 1.5<jupyter_output><empty_output><jupyter_text>Combine the two distributions and define grid of wealth levels associated with $g(w)$<jupyter_code>g = np.append(g_bottom[:-1], g_top)
w = np.append(w_bottom[:-1], tail_vals)<jupyter_output><empty_output><jupyter_text>Plot the distribution of wealth<jupyter_code># Standard 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(w, g, color = '0.4', zorder = 2, label = r'$g(a)$')

ax.set_xlim(-0.3,10)
ax.set_ylim(0,0.35)

ax.set_xlabel(r'Wealth, $w$')
ax.set_ylabel(r'Density, $g(w)$')

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

plt.savefig('../../tex/figs/model.pdf')
# Log-log 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(w, g, color = '0.4', zorder = 2, label = r'$g(a)$')
ax.set(xscale='log', yscale = 'log')

ax.set_xlim(1e-2,1e+3)
ax.set_ylim(1e-8,1)

ax.set_xlabel(r'Wealth, $w$')
ax.set_ylabel(r'Density, $g(w)$')

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

plt.savefig('../../tex/figs/model_log.pdf')<jupyter_output><empty_output><jupyter_text>Simulations Draw from the asset distribution using rejection sampling.First generate the envelope density, $\tilde{g}(w)$, then find $c = \max \frac{g(w)}{\tilde{g}(w)} + \varepsilon$ for some small $\varepsilon$.<jupyter_code># Generate envelope density
density_pareto = funcs.den_pareto2(w, 1.5, np.min(w))

# Find the ratio of the densities over the domain
c = g/density_pareto

# Find the maximal value and add small epsilon
c_max = np.max(c) + 1e-5
print(c_max)<jupyter_output>1.04383186071944<jupyter_text>To draw random numbers from these distributions i need to interpolate between values.I use a linear interpolater for simplicity.<jupyter_code>f = interpolate.interp1d(w,g, bounds_error = False, fill_value = 0)
h = interpolate.interp1d(w,density_pareto, bounds_error = False, fill_value = 0)<jupyter_output><empty_output><jupyter_text>Run the rejection sampling algorithm just to test if everything works.<jupyter_code># Test function
N = 50_000

sample = funcs.act_rejct(N,c_max,h,f)
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(sample, bins = np.linspace(np.min(w),2,500), density = True)
ax.plot(w,g)

ax.set_xlim(np.min(w),2)

plt.show()<jupyter_output><empty_output><jupyter_text>Looks fine **Run Monte Carlo experiment**<jupyter_code># Parameters
N = 10_000
reps = 10_000

# Initialize
G_new = np.empty(shape=reps) # Improved estimator
G_np = np.empty(shape=reps) # Non-parametric estimator

# Loop over reps
for j in range(reps):
    
    # Draw sample using accept-reject
    sample = funcs.act_rejct(N,c_max,h,f)
    
    ## Compute nonparametric gini ## 
    
    G_np[j] = funcs.gini(sample)
    
    ## Compute improved gini ##
    
    # The (1-p)'th percentile
    P = 0.0025
    u = np.percentile(sample,100*(1-P))
    
    # Compute s
    x_low = sample[sample<u] # Is this needed?
    E_low = np.mean(sample[sample<u])
    S = alpha/(alpha-1)*u*P / (alpha/(alpha-1)*u*P + E_low * (1-P))

    # Compute the Gini for top and bottom of sample
    G_high = 1/(2*alpha-1)
    G_low = funcs.gini(sample[sample<u])
    
    # Collect these
    G_new[j] = G_high*P*S + G_low*(1-P)*(1-S) + S - P<jupyter_output><empty_output><jupyter_text>Plot the resulting distribution of Gini coefficient estimators<jupyter_code>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_new, density = True, bins = 150, color = 'red', alpha = 0.7, zorder = 2, label='Semiparametric',  histtype='stepfilled')
ax.hist(G_np, density = True, bins = 150, color = 'blue', alpha = 0.7, zorder = 2, label='Nonparametric',  histtype='stepfilled')

# ax.set_ylim(0,100)
# ax.set_xlim(0.65,1)

ax.set_xlabel(r'Estimated Gini coefficient')
ax.set_ylabel(r'Density')

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

plt.savefig('../../tex/figs/model_estimator.pdf')<jupyter_output><empty_output>