<jupyter_start><jupyter_text>Stationary Equilibrium In this notebook, I illustrate the solution to the **stationary equilibrium** for the housing model.  Setup and Solve<jupyter_code>%load_ext autoreload
%autoreload 2

import time
import numpy as np
import numba as nb
import matplotlib.pyplot as plt
from scipy import optimize
from consav import linear_interp
from scipy import stats

# load the DurableConsumptionModel module
from HousingModel import HousingModelClass

# Set the number of threads in numba
nb.set_num_threads(4)
# Something for figure layout here
plt.rc('font', family='serif')
plt.rc('text', usetex=True)
plt.rc('axes', axisbelow=True)
# Define the model
model = HousingModelClass(name='example_negm',par={'solmethod':'negm', 'do_print':False})

# Solve and simulate
model.precompile_numba()
model.solve()
model.simulate()
model.solve_path()
model.simulate_path()<jupyter_output><empty_output><jupyter_text>Determine aggregate supply so that the equilibrium price is 1<jupyter_code>Hs = np.mean(model.sim.n)
print(Hs)<jupyter_output>0.31576333475798035<jupyter_text>Find stationary equilibrium Find demand curve over a grid of house prices<jupyter_code># Grid over house prices
Nr = 15
p_vec = np.linspace(.2,3.0,Nr)
p_min = p_vec[0]
p_max = p_vec[Nr-1]

# Supply
# Hs = 4.1 # Find some value

# Initialize demand function
Hd = np.zeros(shape=len(p_vec))

# Loop
t0 = time.time()
for i_p,p in enumerate(p_vec):
    
    # Solve and simulate model
    model.par.ph = p
    model.solve()
    model.simulate()
    
    # Aggregate demand
    Hd[i_p] = np.mean(model.sim.d[175:,:])
    
    # Bound on p
    diff = Hs - Hd[i_p]
    if diff < 0: p_min = np.fmax(p_min, p_vec[i_p])
    if diff > 0: p_max = np.fmin(p_max, p_vec[i_p])
    
    # print(i_p)
    
t1 = time.time()
print(t1-t0)
fig = plt.figure(figsize=(6,4))
ax = fig.add_subplot(1,1,1)

ax.plot(Hd,p_vec, color = 'red', label = 'Demand',zorder=2)
ax.scatter(Hd,p_vec, color = 'red',zorder=3, marker='o', s=10)

ax.vlines(0.32,0,10, color = 'blue', linestyle = '-', label ='Supply')
# ax.plot(Hd,supply, color = 'blue', linestyle = ':')

ax.hlines(p_min,0,2, linestyle='--', color = '0.2', linewidth = 1)
ax.hlines(p_max,0.0,2, linestyle='--', color = '0.2', linewidth = 1)

ax.set_xlim(0,.68)
ax.set_ylim(.1,3.1)
ax.set_xlabel('Quantity')
ax.set_ylabel('Price')

ax.grid(visible=True, which = 'major', linestyle='-', linewidth=0.5, color='0.9')
ax.tick_params(axis='both',bottom=True,top=True,left=True,right=True,direction='in',which='both')
ax.legend(frameon = True, edgecolor = 'k', facecolor = 'white', framealpha=1, fancybox=False)

plt.savefig('../../tex/figs/equilibrium.pdf')<jupyter_output><empty_output><jupyter_text>Last step<jupyter_code>def obj(p,model):
    
    model.par.ph = p
    model.solve()
    model.simulate()
    return np.mean(model.sim.n[-1,:]) - Hs

t0 = time.time()

opt = optimize.root_scalar(obj,bracket=[p_min,p_max],method='bisect',args=(model,),maxiter=100)
model.par.ph = opt.root
assert opt.converged

t1 = time.time()
print(t1-t0)
print(opt)<jupyter_output><empty_output><jupyter_text>Plot Results Plot the **policy functions** and **stationary distributions** from the model.Also illustrate convergence of the distributions  Plot the policy functions Non-durable Consumption<jupyter_code># Initialize the figure
fig = plt.figure(figsize=(4,3))
ax = fig.add_subplot(1,1,1)

# Plot the avalanche sizes
ax.plot(model.par.grid_m, model.sol.c_keep[0,0,0,9,:], color = '0', linewidth = '1.5', label = r'$y=y_1$')
ax.plot(model.par.grid_m, model.sol.c_keep[0,0,1,9,:], color = '0.2', linewidth = '1.5', label = r'$y=y_2$')
ax.plot(model.par.grid_m, model.sol.c_keep[0,0,2,9,:], color = '0.4', linewidth = '1.5', label = r'$y=y_3$')
ax.plot(model.par.grid_m, model.sol.c_keep[0,0,3,9,:], color = '0.6', linewidth = '1.5', label = r'$y=y_4$')
ax.plot(model.par.grid_m, model.sol.c_keep[0,0,4,9,:], color = '0.8', linewidth = '1.5', label = r'$y=y_5$')

# Plot setting
ax.set_xlabel('Cash on Hand')
ax.set_ylabel('Consumption')

ax.set_xlim(model.par.grid_m[0],model.par.grid_m[-1])
ax.set_ylim(0,1.7)

ax.grid(visible=True, which = 'major', linestyle='-', linewidth=0.5, color='0.9')
ax.tick_params(axis='both',bottom=True,top=True,left=True,right=True,direction='in',which='both')
ax.legend(frameon = True, edgecolor = 'k', facecolor = 'white', framealpha=1, fancybox=False)

plt.savefig('../../tex/figs/c_keep.pdf')
# Initialize the figure
fig = plt.figure(figsize=(4,3))
ax = fig.add_subplot(1,1,1)

# Plot the avalanche sizes
ax.plot(model.par.grid_x, model.sol.c_adj[0,0,0,:], linestyle = '-', color = '0', linewidth = '1.5', label = r'$y=y_1$')
ax.plot(model.par.grid_x, model.sol.c_adj[0,0,1,:], linestyle = '-', color = '0.2', linewidth = '1.5', label = r'$y=y_2$')
ax.plot(model.par.grid_x, model.sol.c_adj[0,0,2,:], linestyle = '-', color = '0.4', linewidth = '1.5', label = r'$y=y_3$')
ax.plot(model.par.grid_x, model.sol.c_adj[0,0,3,:], linestyle = '-', color = '0.6', linewidth = '1.5', label = r'$y=y_4$')
ax.plot(model.par.grid_x, model.sol.c_adj[0,0,4,:], linestyle = '-', color = '0.8', linewidth = '1.5', label = r'$y=y_5$')

# Plot setting
ax.set_xlabel('Cash on Hand')
ax.set_ylabel('Consumption')

ax.set_xlim(0,15)
# ax.set_ylim(0,1.35)

ax.grid(visible=True, which = 'major', linestyle='-', linewidth=0.5, color='0.9')
ax.tick_params(axis='both',bottom=True,top=True,left=True,right=True,direction='in',which='both')
ax.legend(frameon = True, edgecolor = 'k', facecolor = 'white', framealpha=1, fancybox=False)

plt.savefig('../../tex/figs/c_adj.pdf')<jupyter_output><empty_output><jupyter_text>Housing demand<jupyter_code># Initialize the figure
fig = plt.figure(figsize=(4,3))
ax = fig.add_subplot(1,1,1)

# Plot the avalanche sizes
ax.plot(model.par.grid_x, model.sol.d_adj[0,0,0,:], linestyle = '-', color = '0', linewidth = '1.5', label = r'$y=y_1$')
ax.plot(model.par.grid_x, model.sol.d_adj[0,0,1,:], linestyle = '-', color = '0.2', linewidth = '1.5', label = r'$y=y_2$')
ax.plot(model.par.grid_x, model.sol.d_adj[0,0,2,:], linestyle = '-', color = '0.4', linewidth = '1.5', label = r'$y=y_3$')
ax.plot(model.par.grid_x, model.sol.d_adj[0,0,3,:], linestyle = '-', color = '0.6', linewidth = '1.5', label = r'$y=y_4$')
ax.plot(model.par.grid_x, model.sol.d_adj[0,0,4,:], linestyle = '-', color = '0.8', linewidth = '1.5', label = r'$y=y_5$')

# Plot setting
ax.set_xlabel('Cash on Hand')
ax.set_ylabel('Housing Demand')

ax.set_xlim(0,20)
ax.set_ylim(0,1.5)

ax.grid(visible=True, which = 'major', linestyle='-', linewidth=0.5, color='0.9')
ax.tick_params(axis='both',bottom=True,top=True,left=True,right=True,direction='in',which='both')
ax.legend(frameon = True, edgecolor = 'k', facecolor = 'white', framealpha=1, fancybox=False)

plt.savefig('../../tex/figs/h_adj.pdf')<jupyter_output><empty_output><jupyter_text>Convergence of dist<jupyter_code># plt.plot(range(200,0,-1),model1.sol.dist)

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

# Plot and scatter
ax.plot(range(100,0,-1),model.sol.dist[100:], color = 'blue', linewidth = '1.5')

# Plot setting
ax.set_xlabel('Iteration Number')
ax.set_ylabel('Distance')

ax.set_xlim(0,100)
# ax.set_ylim(0.25,2.0)

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

plt.savefig('../../tex/figs/dist.pdf')<jupyter_output><empty_output><jupyter_text>Plot Distributions  Wealth distribution First, plot the distribution in a normal plot<jupyter_code># Define wealth
w = (model.sim.m + model.par.ph*model.sim.n) / model.sim.y
# w2 = model.sim.m 

# w = model.sim.a + model.sim.d

# Initialize the figure
fig = plt.figure(figsize=(4,3))
ax = fig.add_subplot(1,1,1)

# Plot the avalanche sizes
# ax.hist(w[-1,:], bins = 100, density = True, edgecolor = 'k', linewidth=.5, zorder=2, color ='royalblue', alpha = 1)
ax.hist(w[-1,:], bins = np.linspace(0,25,50), density = True, edgecolor = 'k', linewidth=.5, zorder=2, color ='royalblue', alpha = 1)

# ax.hist(w2[-1,:], bins = 50, density = True, edgecolor = 'k', linewidth=.5, zorder=2, color ='royalblue', alpha = .5)

# Plot setting
ax.set_xlabel('Relative Wealth')
ax.set_ylabel('Density')
ax.set_xlim(0,25)

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

plt.savefig('../../tex/figs/w_distr.pdf')<jupyter_output><empty_output><jupyter_text>Use a log-log plot<jupyter_code>## Define wealth
w = (model.sim.m + model.par.ph*model.sim.n) / model.sim.y

# Initialize the figure
fig = plt.figure(figsize=(4,3))
ax = fig.add_subplot(1,1,1)

# Plot the avalanche sizes
ax.hist(w[-1,:], np.logspace(.1,2,60), density = True, edgecolor = 'k', linewidth=.5, zorder=2, color ='royalblue')

# Plot setting
ax.set_xlabel('Relative Wealth')
ax.set_ylabel('Density')

ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlim(1,80)

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

plt.savefig('../../tex/figs/w_distr_log.pdf')<jupyter_output><empty_output><jupyter_text>Aggregate housing demand<jupyter_code># Initialize the figure
fig = plt.figure(figsize=(4,3))
ax = fig.add_subplot(1,1,1)

# Plot the avalanche sizes
ax.hist(model.sim.d[-1,:], model.par.grid_n, density = False, facecolor='orangered', edgecolor = 'k', zorder=2)

# Plot setting
ax.set_xlabel('Housing')
ax.set_ylabel('Nr. of Households')

ax.set_xlim(0,1.5)

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

plt.savefig('../../tex/figs/h_distr.pdf')<jupyter_output><empty_output><jupyter_text>Plot Convergence of Simulation  Simulate and check for convergence<jupyter_code># Vector of aggregates
h_agg = np.mean(model.sim.d, axis=1)
n_agg = np.mean(model.sim.n, axis=1)
c_agg = np.mean(model.sim.c, axis=1)
a_agg = np.mean(model.sim.a, axis=1)
m_agg = np.mean(model.sim.m, axis=1)
y_agg = np.mean(model.sim.y, axis=1)<jupyter_output><empty_output><jupyter_text>Divide time series with the last period value to illustrate convergence<jupyter_code>h_conv = h_agg / h_agg[-1]
c_conv = c_agg / c_agg[-1]
a_conv = a_agg / a_agg[-1]
# Initialize the figure
fig = plt.figure(figsize=(6,4))
ax = fig.add_subplot(1,1,1)

# Plot the avalanche sizes
ax.plot(range(model.par.sim_T),h_conv, linestyle = '-', color = 'purple', linewidth = '1.1', label= 'Housing')
ax.plot(range(model.par.sim_T),c_conv, linestyle = '-', color = 'green', linewidth = '1.1', label='Non-durable')
ax.plot(range(model.par.sim_T),a_conv, linestyle = '-', color = 'blue', linewidth = '1.1', label='End of period assets')

ax.hlines(1,0,200, linestyle ='--', color = '0.2', linewidth = 1.3)

# Plot setting
ax.set_xlabel('Period Number')
ax.set_ylabel('Deviation From Steady State')

ax.set_xlim(0,200)
ax.set_ylim(0.8,1.2)

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

plt.savefig('../../tex/figs/convergence.pdf')<jupyter_output><empty_output><jupyter_text>Descriptive Statistics from SS Compute variables used<jupyter_code># Housing share of total consumption
h_share = model.sim.n[-1,:] / ((model.sim.m[-1,:] + model.sim.n[-1,:]) / model.sim.y[-1,:])
# try binned scatterplot
h_demand_means, bin_edges, binnumber = stats.binned_statistic(model.sim.m[-1,:],model.sim.d[-1,:], bins = np.logspace(0,2,30), statistic = 'mean')
h_share_means = stats.binned_statistic(model.sim.m[-1,:],h_share, bins = np.logspace(0,2,30), statistic = 'mean').statistic

# Compute bins
bin_width = (bin_edges[1] - bin_edges[0])
bin_centers = bin_edges[1:] - bin_width/2<jupyter_output><empty_output><jupyter_text>Binscatter of housing demand as a function of beginning of period cash on hand<jupyter_code># Initialize the figure
fig = plt.figure(figsize=(6,4))
ax = fig.add_subplot(1,1,1)

# Plot the avalanche sizes
ax.scatter(bin_centers, h_demand_means, s=50, c='orange', marker='o', edgecolors = '0.2', zorder=2) 

# Plot setting
ax.set_xlabel('Cash-On-Hand')
ax.set_ylabel('Demand')

ax.set_xlim(0,17.5)
# ax.set_ylim(0,1.4)

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

plt.savefig('../../tex/figs/h_demand_bins.pdf')<jupyter_output><empty_output><jupyter_text>Binscatter of the share of housing in total consumption as a function of beginning of period cash on hand<jupyter_code># Initialize the figure
fig = plt.figure(figsize=(6,4))
ax = fig.add_subplot(1,1,1)

# Plot the avalanche sizes
ax.scatter(bin_centers, h_share_means, s=50, c='magenta', marker='o', edgecolors = '0.2', zorder=2) 

# Plot setting
ax.set_xlabel('Cash-On-Hand')
ax.set_ylabel('Share')

ax.set_xlim(0,16.4)
# ax.set_ylim(0.07,0.3)

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

plt.savefig('../../tex/figs/h_share_bins.pdf')<jupyter_output><empty_output><jupyter_text>Find SS for lower interest rate Save some statistics needed<jupyter_code>h_init = model.sim.n
w_init = (model.sim.m + model.par.ph*model.sim.n) / model.sim.y<jupyter_output><empty_output><jupyter_text>Define model2 with lower interest rates<jupyter_code># Define the model
R_low = 1.02
model2 = HousingModelClass(name='example_negm',par={'solmethod':'negm', 'R':R_low, 'ph':p, 'do_print':False})

# Solve and simulate
model2.precompile_numba()
model2.solve()
model2.simulate()
print(Hs)
# Grid over house prices
Nr = 15
p_vec = np.linspace(.3,3.0,Nr)
p_min = p_vec[0]
p_max = p_vec[Nr-1]

# Supply set to same value as the previous model

# Initialize demand function
Hd = np.zeros(shape=len(p_vec))

# Loop
t0 = time.time()
for i_p,p in enumerate(p_vec):
    
    # Solve and simulate model
    model2.par.ph = p
    model2.solve()
    model2.simulate()
    
    # Aggregate demand
    Hd[i_p] = np.mean(model2.sim.d[-1,:])
    
    # Bound on p
    diff = Hs - Hd[i_p]
    if diff < 0: p_min = np.fmax(p_min, p_vec[i_p])
    if diff > 0: p_max = np.fmin(p_max, p_vec[i_p])
    
    # print(i_p)
    
t1 = time.time()
def obj(p,model2):
    
    model2.par.ph = p
    model2.solve()
    model2.simulate()
    return np.mean(model2.sim.d[-1,:]) - Hs

t0 = time.time()

opt = optimize.root_scalar(obj,bracket=[p_min,p_max],method='bisect',args=(model2,),maxiter=100)
model2.par.ph = opt.root
assert opt.converged

t1 = time.time()
print(t1-t0)
opt.root
# Define wealth
w = (model2.sim.m + model2.par.ph*model2.sim.n) / model2.sim.y

# w = model.sim.a + model.sim.d

# Initialize the figure
fig = plt.figure(figsize=(4,3))
ax = fig.add_subplot(1,1,1)

# Plot the avalanche sizes
ax.hist(w[-1,:], bins = np.linspace(0,25,50), density = True, edgecolor = 'k', linewidth=.5, zorder=2, color ='royalblue', alpha = .5, label = r'Low $R$')
ax.hist(w_init[-1,:], bins = np.linspace(0,25,50), density = True, edgecolor = 'k', linewidth=.5, zorder=2, color ='green', alpha = .5, label = r'High $R$')

# Plot setting
ax.set_xlabel('Relative Wealth')
ax.set_ylabel('Density')
ax.set_xlim(0,25)

ax.grid(visible=True, which = 'major', linestyle='-', linewidth=0.5, color='0.9')
ax.tick_params(axis='both',bottom=True,top=True,left=True,right=True,direction='in',which='both')
ax.legend(frameon = True, edgecolor = 'k', facecolor = 'white', framealpha=1, fancybox=False, prop={'size': 10})

plt.savefig('../../tex/figs/w_distr_compare.pdf')
# Initialize the figure
fig = plt.figure(figsize=(4,3))
ax = fig.add_subplot(1,1,1)

# Plot the avalanche sizes
ax.hist(model2.sim.d[-1,:], model2.par.grid_n, density = False, facecolor='royalblue', edgecolor = 'k', zorder=2, alpha = .4, label = r'Low $R$')
ax.hist(h_init[-1,:], model2.par.grid_n, density = False, facecolor='green', edgecolor = 'k', zorder=2, alpha = .4, label = r'High $R$')

# Plot setting
ax.set_xlabel('Housing')
ax.set_ylabel('Probability')

ax.set_xlim(0,1.5)

ax.grid(visible=True, which = 'major', linestyle='-', linewidth=0.5, color='0.9')
ax.tick_params(axis='both',bottom=True,top=True,left=True,right=True,direction='in',which='both')
ax.legend(frameon = True, edgecolor = 'k', facecolor = 'white', framealpha=1, fancybox=False,  prop={'size': 10})

plt.savefig('../../tex/figs/h_distr_compare.pdf')
# try binned scatterplot
h_demand_means, bin_edges, binnumber = stats.binned_statistic(model2.sim.m[-1,:],model2.sim.d[-1,:], bins = np.logspace(0,2,30), statistic = 'mean')

# Compute bins
bin_width = (bin_edges[1] - bin_edges[0])
bin_centers = bin_edges[1:] - bin_width/2
# Initialize the figure
fig = plt.figure(figsize=(6,4))
ax = fig.add_subplot(1,1,1)

# Plot the avalanche sizes
ax.scatter(bin_centers, h_demand_means, s=50, c='red', marker='o', edgecolors = '0.2', zorder=2) 

# Plot setting
ax.set_xlabel('Cash-On-Hand')
ax.set_ylabel('Demand')

ax.set_xlim(0,17.5)
ax.set_ylim(0,1.2)

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

plt.savefig('../../tex/figs/h_demand_bins_low.pdf')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19431198910081743, 'fraction_numerical': 0.03371934604904632, 'mean_word_length': 5.59812734082397, 'pattern_counts': {'":': 0, '<': 70, '<?xml version=': 0, '>': 70, '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': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '373926', 'license': 'MIT', 'n_tokens_mistral': 7467, 'n_tokens_neox': 7267, 'n_words': 1273, 'path': '00 - Stationary.ipynb', 'repo_name': 'madsankern/MA_thesis', 'size': '17616'}