In [1]:
%matplotlib inline
In [2]:
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sn
sn.set_style("whitegrid")
import sympy as sym
import pycollocation
import pyam
In [3]:
pyam.__version__
Out[3]:
In [4]:
# define some workers skill
x, loc1, mu1, sigma1 = sym.var('x, loc1, mu1, sigma1')
skill_cdf = 0.5 + 0.5 * sym.erf((sym.log(x - loc1) - mu1) / sym.sqrt(2 * sigma1**2))
skill_params = {'loc1': 1e0, 'mu1': 0.0, 'sigma1': 1.0}
workers = pyam.Input(var=x,
cdf=skill_cdf,
params=skill_params,
bounds=(1.2, 1e1), # guesses for the alpha and (1 - alpha) quantiles!
alpha=0.005,
measure=1.0
)
# define some firms
y, loc2, mu2, sigma2 = sym.var('y, loc2, mu2, sigma2')
productivity_cdf = 0.5 + 0.5 * sym.erf((sym.log(y - loc2) - mu2) / sym.sqrt(2 * sigma2**2))
productivity_params = {'loc2': 1e0, 'mu2': 0.0, 'sigma2': 1.0}
firms = pyam.Input(var=y,
cdf=productivity_cdf,
params=productivity_params,
bounds=(1.2, 1e1), # guesses for the alpha and (1 - alpha) quantiles!
alpha=0.005,
measure=1.0
)
Note that we are shifting the distributions of worker skill and firm productivity to the right by 1.0 in order to try and avoid issues with having workers (firms) with near zero skill (productivity).
In [5]:
xs = np.linspace(workers.lower, workers.upper, 1e4)
plt.plot(xs, workers.evaluate_pdf(xs))
plt.xlabel('Worker skill, $x$', fontsize=20)
plt.show()
In [6]:
# define symbolic expression for CES between x and y
omega_A, sigma_A = sym.var('omega_A, sigma_A')
A = ((omega_A * x**((sigma_A - 1) / sigma_A) +
(1 - omega_A) * y**((sigma_A - 1) / sigma_A))**(sigma_A / (sigma_A - 1)))
# define symbolic expression for CES between x and y
r, l, omega_B, sigma_B = sym.var('r, l, omega_B, sigma_B')
B = ((omega_B * r**((sigma_B - 1) / sigma_B) +
(1 - omega_B) * l**((sigma_B - 1) / sigma_B))**(sigma_B / (sigma_B - 1)))
F = A * B
In [7]:
# negative assortativity requires that sigma_A * sigma_B > 1
F_params = {'omega_A':0.25, 'omega_B':0.5, 'sigma_A':2.0, 'sigma_B':1.0 }
In [8]:
problem = pyam.AssortativeMatchingProblem(assortativity='negative',
input1=workers,
input2=firms,
F=sym.limit(F, sigma_B, 1),
F_params=F_params)
In [9]:
solver = pycollocation.OrthogonalPolynomialSolver(problem)
Currently I guess that $\mu(x)$ is has the form...
$$ \hat{\mu}(x) = \beta_0 + \beta_1 f(x) $$(i.e., a linear translation) of some function $f$. Using my $\hat{\mu}(x)$, I can then back out a guess for $\theta(x)$ implied by the model...
$$ \hat{\theta}(x) = \frac{H(x)}{\hat{\mu}'(x)} $$
In [59]:
initial_guess = pyam.OrthogonalPolynomialInitialGuess(solver)
initial_polys = initial_guess.compute_initial_guess("Chebyshev",
degrees={'mu': 40, 'theta': 70},
f=lambda x, alpha: x**alpha,
alpha=1.0)
In [60]:
# quickly plot the initial conditions
xs = np.linspace(workers.lower, workers.upper, 1000)
plt.plot(xs, initial_polys['mu'](xs))
plt.plot(xs, initial_polys['theta'](xs))
plt.grid('on')
In [61]:
domain = [workers.lower, workers.upper]
initial_coefs = {'mu': initial_polys['mu'].coef,
'theta': initial_polys['theta'].coef}
solver.solve(kind="Chebyshev",
coefs_dict=initial_coefs,
domain=domain,
method='hybr')
In [62]:
solver.result.success
Out[62]:
In [63]:
viz = pyam.Visualizer(solver)
In [64]:
viz.interpolation_knots = np.linspace(workers.lower, workers.upper, 1000)
viz.residuals.plot()
plt.show()
In [58]:
viz.normalized_residuals[['mu', 'theta']].plot(logy=True)
plt.show()
In [53]:
viz.solution.tail()
Out[53]:
In [25]:
viz.solution[['mu', 'theta']].plot(subplots=True)
plt.show()
In [26]:
viz.solution[['Fxy', 'Fyl']].plot()
plt.show()
In [27]:
viz.solution[['factor_payment_1', 'factor_payment_2']].plot(subplots=True)
plt.show()
In [28]:
fig, axes = plt.subplots(1, 2, sharey=True)
axes[0].scatter(viz.solution.factor_payment_1, viz.solution.theta, alpha=0.5,
edgecolor='none')
axes[0].set_ylim(0, 1.05 * viz.solution.theta.max())
axes[0].set_xlabel('Wages, $w$')
axes[0].set_ylabel(r'Firm size, $\theta$')
axes[1].scatter(viz.solution.factor_payment_2, viz.solution.theta, alpha=0.5,
edgecolor='none')
axes[1].set_xlabel(r'Profits, $\pi$')
plt.show()
In [29]:
# to get correlation just use pandas!
viz.solution.corr()
Out[29]:
In [30]:
# or a subset
viz.solution[['theta', 'factor_payment_1']].corr()
Out[30]:
In [31]:
# or actual values!
viz.solution.corr().loc['theta']['factor_payment_1']
Out[31]:
As you can see, the theta function is hump-shaped. Nothing special, but when calculating the pdf some arrangements have to be done for this: sort the thetas preserving the order (so we can relate them to their xs) and then use carefully the right x for calculating the pdf.
The principle of Philipp's trick is:
$pdf_x(x_i)$ can be interpreted as number of workers with ability x. $\theta_i$ is the size of the firms that employs workers of kind $x_i$. As all firms that match with workers type $x_i$ choose the same firm size, $pdf_x(x_i)/\theta_i$ is the number of firms of size $\theta_i$.
Say there are 100 workers with ability $x_i$, and their associated firm size $\theta_i$ is 2. Then there are $100/2 = 50$ $ \theta_i$ firms
In [32]:
fig, axes = plt.subplots(1, 3)
theta_pdf = viz.compute_pdf('theta', normalize=True)
theta_pdf.plot(ax=axes[0])
axes[0].set_xlabel(r'Firm size, $\theta$')
axes[0].set_title(r'pdf')
theta_cdf = viz.compute_cdf(theta_pdf)
theta_cdf.plot(ax=axes[1])
axes[1].set_title(r'cdf')
axes[1].set_xlabel(r'Firm size, $\theta$')
theta_sf = viz.compute_sf(theta_cdf)
theta_sf.plot(ax=axes[2])
axes[2].set_title(r'sf')
axes[2].set_xlabel(r'Firm size, $\theta$')
plt.tight_layout()
plt.show()
In [33]:
fig, axes = plt.subplots(1, 3)
factor_payment_1_pdf = viz.compute_pdf('factor_payment_1', normalize=True)
factor_payment_1_pdf.plot(ax=axes[0])
axes[0].set_title(r'pdf')
factor_payment_1_cdf = viz.compute_cdf(factor_payment_1_pdf)
factor_payment_1_cdf.plot(ax=axes[1])
axes[1].set_title(r'cdf')
factor_payment_1_sf = viz.compute_sf(factor_payment_1_cdf)
factor_payment_1_sf.plot(ax=axes[2])
axes[2].set_title(r'sf')
plt.tight_layout()
plt.show()
In [34]:
fig, axes = plt.subplots(1, 3)
factor_payment_2_pdf = viz.compute_pdf('factor_payment_2', normalize=True)
factor_payment_2_pdf.plot(ax=axes[0])
axes[0].set_title(r'pdf')
factor_payment_2_cdf = viz.compute_cdf(factor_payment_2_pdf)
factor_payment_2_cdf.plot(ax=axes[1])
axes[1].set_title(r'cdf')
factor_payment_2_sf = viz.compute_sf(factor_payment_2_cdf)
factor_payment_2_sf.plot(ax=axes[2])
axes[2].set_title(r'sf')
plt.tight_layout()
plt.show()
In [54]:
from IPython.html import widgets
In [55]:
def interactive_plot(viz, omega_A=0.25, omega_B=0.5, sigma_A=0.5, sigma_B=1.0,
loc1=1.0, mu1=0.0, sigma1=1.0, loc2=1.0, mu2=0.0, sigma2=1.0):
# update new parameters as needed
new_F_params = {'omega_A': omega_A, 'omega_B': omega_B,
'sigma_A': sigma_A, 'sigma_B': sigma_B}
viz.solver.problem.F_params = new_F_params
new_input1_params = {'loc1': loc1, 'mu1': mu1, 'sigma1': sigma1}
viz.solver.problem.input1.params = new_input1_params
new_input2_params = {'loc2': loc2, 'mu2': mu2, 'sigma2': sigma2}
viz.solver.problem.input2.params = new_input2_params
# solve the model using a hotstart initial guess
domain = [viz.solver.problem.input1.lower, viz.solver.problem.input1.upper]
initial_coefs = viz.solver._coefs_array_to_dict(viz.solver.result.x, viz.solver.degrees)
viz.solver.solve(kind="Chebyshev",
coefs_dict=initial_coefs,
domain=domain,
method='hybr')
if viz.solver.result.success:
viz._Visualizer__solution = None # should not need to access this!
viz.interpolation_knots = np.linspace(domain[0], domain[1], 1000)
viz.solution[['mu', 'theta']].plot(subplots=True)
viz.normalized_residuals[['mu', 'theta']].plot(logy=True)
else:
print "Foobar!"
In [56]:
viz_widget = widgets.fixed(viz)
# widgets for the model parameters
eps = 1e-2
omega_A_widget = widgets.FloatSlider(value=0.25, min=eps, max=1-eps, step=eps,
description=r"$\omega_A$")
sigma_A_widget = widgets.FloatSlider(value=0.5, min=eps, max=1-eps, step=eps,
description=r"$\sigma_A$")
omega_B_widget = widgets.FloatSlider(value=0.5, min=eps, max=1-eps, step=eps,
description=r"$\omega_B$")
sigma_B_widget = widgets.fixed(1.0)
# widgets for input distributions
loc_widget = widgets.fixed(1.0)
mu_1_widget = widgets.FloatSlider(value=0.0, min=-1.0, max=1.0, step=eps,
description=r"$\mu_1$")
mu_2_widget = widgets.FloatSlider(value=0.0, min=-1.0, max=1.0, step=eps,
description=r"$\mu_2$")
sigma_1_widget = widgets.FloatSlider(value=1.0, min=eps, max=2-eps, step=eps,
description=r"$\sigma_1$")
sigma_2_widget = widgets.FloatSlider(value=1.0, min=eps, max=2-eps, step=eps,
description=r"$\sigma_2$")
In [57]:
widgets.interact(interactive_plot, viz=viz_widget, omega_A=omega_A_widget,
sigma_A=sigma_A_widget, omega_B=omega_B_widget,
sigma_B=sigma_B_widget, sigma1=sigma_1_widget,
loc1=loc_widget, mu1 = mu_1_widget,
loc2=loc_widget, sigma2=sigma_2_widget, mu2 = mu_2_widget)
In [34]:
# widget is changing the parameters of the underlying solver
solver.result.x
Out[34]:
In [ ]: