# Import packages
import numpy as np # Required for vectorized operations, linear algebra, random number generation and probability distributions
import pandas as pd # Required for data handling
import scipy.stats # Required for statistical tests, functions and distributions
import matplotlib.pyplot as plt # Required for plotting results
# Set parameters (N = 20,000; M = 20,000)
N = 20000 # Number of workers
M = 20000 # Number of firms
r = 0.02 # Interest rate
gamma = 0.05 # Rate of workerfirm matching (1/gamma is average duration of job)
beta = 0.05 # Rate of worker search effort adjustment (1/beta is average period length over time (size 100)
E_array = np.empty(100) # Array to store aggregate employment over time (size 100)
U_array = np.empty(100) # Array to store aggregate unemployment over time (size 100)
UR_array = np.empty(100) # Array to store aggregate unemployment rate over time (size 100)
P_array = np.empty(100) # Array to store aggregate search effort of workers over time (size 100)
Q_array = np.empty(100) # Array to store aggregate search effort of firms over time (size 100)
MPW_array = np.empty(100) # Array to store aggregate marginal product of workers over time (size 100)
MPF_array = np.empty(100) # Array to store aggregate marginal product of firms over time (size 100)
VR_array = np.empty(100) # Array to store aggregate vacancy rate over time (size 100)
sigma = 1 # Standard deviation of productivity shocks
alpha = 0.15 # Rate of firm search effort adjustment (1/alpha is average period length over time (size 100))
# Generate worker and firm productivity shocks
theta = np.random.normal(0, 2, N) # Worker productivity shock; Normal distribution with mean 0 and standard deviation 2
psi = np.random.normal(0, 2, M) # Firm productivity shock; Normal distribution with mean 0 and standard deviation 2
lambda_ij = np.exp((theta[:,np.newaxis]psi[np.newaxis,:])**2/(2*sigma**2)) # Matching productivity
Lambda_ij = lambda_ij/(1+lambda_ij) # Probability of workerfirm match
delta = 1/(1+lambda_ij) # Probability of workerfirm mismatch
# Initialize equilibrium aggregate search effort values (bar)
P = N/M*Lambda_ij.sum(axis=1).mean() # Equilibrium aggregate search effort for workers (average number of job applications sent out by all workers per period)
Q = M/N*Lambda_ij.sum(axis=0).mean() # Equilibrium aggregate search effort for firms (average number of vacancies posted by all firms per period)
# Initialize search efforts based on productivity shocks and equilibrium aggregate search effort values (bar)
S = P*Lambda_ij/Lambda_ij.sum(axis=1)[:,np.newaxis] # Worker search effort; Number of job applications sent out by each worker per period
s = Q*Lambda_ij/Lambda_ij.sum(axis=0)[np.newaxis,:] # Firm search effort; Number of vacancies posted by each firm per period
# Calculate equilibrium aggregate search effort values (bar)
P = S.sum() # Equilibrium aggregate search effort for workers (average number of job applications sent out by all workers per period)
Q = s.sum() # Equilibrium aggregate search effort for firms (average number of vacancies posted by all firms per period)
# Calculate aggregate marginal products, vacancy rates, employment and unemployment using equilibrium aggregate search effort values (bar)
MPW = 1S/P # Aggregate marginal product of workers (average wages paid to all workers per period)
MPF = 1s/Q # Aggregate marginal product of firms (average profits generated by all firms per period)
VR = s/Q # Aggregate vacancy rate (average number of vacancies posted by all firms per period over average number of vacancies plus average number of jobs created by all firms per period)
E = MPW/(1+r)*gamma/(1gamma) # Aggregate employment (average number of jobs created by all firms per period)
U = MPW/r*(1gamma)*(1MPW/r)*gamma/(1gamma) # Aggregate unemployment (average number of job applications received by all firms per period)
UR = U/(E+U) # Aggregate unemployment rate (average number of job applications received by all firms per period over average number of vacancies plus average number of jobs created by all firms per period)
# Store results in arrays
E_array[0] = E
U_array[0] = U
UR_array[0] = UR
P_array[0] = P
Q_array[0] = Q
MPW_array[0] = MPW
MPF_array[0] = MPF
VR_array[0] = VR
# Simulate model for 100 periods (t)
for t in range(1,100):
# Calculate change in worker and firm search efforts based on previous period's equilibrium aggregate search effort values (bar)
dP = beta*((Q/N)P) # Change in worker search effort; Change in average number of job applications sent out by all workers per period
dQ = alpha*((P/M)Q) # Change in firm search effort; Change in average number of vacancies posted by all firms per period
# Update equilibrium aggregate search effort values (bar)
P = P+dP # Equilibrium aggregate search effort for workers (average number of job applications sent out by all workers per period)
Q = Q+dQ # Equilibrium aggregate search effort for firms (average number of vacancies posted by all firms per period)
# Calculate aggregate marginal products, vacancy rates, employment and unemployment using equilibrium aggregate search effort values (bar)
MPW = 1S/P # Aggregate marginal product of workers (average wages paid to all workers per period)
MPF = 1s/Q # Aggregate marginal product of firms (average profits generated by all firms per period)
VR = s/Q # Aggregate vacancy rate (average number of vacancies posted by all firms per period over average number of vacancies plus average number of jobs created by all firms per period)
E = MPW/(1+r)*gamma/(1gamma) # Aggregate employment (average number of jobs created by all firms per period)
U = MPW/r*(1gamma)*(1MPW/r)*gamma/(1gamma) # Aggregate unemployment (average number of job applications received by all firms per period)
UR = U/(E+U) # Aggregate unemployment rate (average number of job applications received by all firms per period over average number of vacancies plus average number of jobs created by all firms per period)
# Store results in arrays
E_array[t] = E
U_array[t] = U
UR_array[t] = UR
P_array[t] = P
Q_array[t] = Q
MPW_array[t] = MPW
MPF_array[t] = MPF
VR_array[t] = VR
# Plot simulation results
fig, axs = plt.subplots(2,2, figsize=(15,7), sharex=True) # Set up subplots in figure with size 15 x 7 inches; Share xaxis between subplots
axs[0,0].plot(E_array) # Plot aggregate employment
axs[0,0].set_title('Aggregate Employment') # Set title for subplot
axs[0,1].plot(U_array) # Plot aggregate unemployment
axs[0,1].set_title('Aggregate Unemployment') # Set title for subplot
axs[1,0].plot(UR_array) # Plot aggregate unemployment rate
axs[1,0].set_title('Aggregate Unemployment Rate') # Set title for subplot
axs[1,1].plot(VR_array) # Plot aggregate vacancy rate
axs[1,1].set_title('Aggregate Vacancy Rate') # Set title for subplot
plt.show() # Display figure
KernelInterrupted: Execution interrupted by the Jupyter kernel.
