I managed to modify the original version just adding new probabilities: φ, σ, γ, β and transforming the original matrices in 3x3 adding EUT as new variable.
However, I have a problem. I am not able to build the key matrice of the mode A_hat such that it has eigenvalues lower than 1, I think for the last part of the code we need a a diagonalizable matrix but I do not know how to do it. I try to modify the value of the various probabilities but it didn’t really changed. Am I doing something wrong? Is that the case I cannot obtain it?
%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (11, 5) #set default figure size
import numpy as np
from quantecon import MarkovChain
from scipy.stats import norm
from scipy.optimize import brentq
from quantecon.distributions import BetaBinomial
from numba import jit
class NewLakeModel:
"""
We compute dynamics of unemployment, employment and employment
under trade union.
parameters are:
λ : scalar
The job finding rate for currently unemployed workers
α : scalar
The dismissal rate for currently employed workers
b : scalar
Entry rate into the labor force
d : scalar
Exit rate from the labor force
φ : scalar
job finding rate for currently unemployed workers toward a job with TU
σ : scalar
The dismissal rate for currently employed workers under TU
γ : scalar
Switch rate for employed workers under TU towards job without TU
β : scalar
Switch rate for employed workers to be employed under a TU
"""
def __init__(self, λ=0.224,α=0.013, b=0.0124, d=0.0082, φ=0.01,σ=0.02,γ=0.03, β=0.081):
self.λ, self.α, self.b, self.d, self.φ, self.σ, self.γ, self.β = λ, α, b, d, φ, σ, γ, β
λ, α, b, d, φ, σ, c, β = self.λ, self.α, self.b, self.d, self.φ, self.σ, self.γ, self.β
self.g = b - d
self.A = np.array([[(1-d)*(1-λ)*(1-φ) + b, (1-d)*(1-β)*α + b, (1-d)*σ*(1-γ) + b],
[(1-d)*(1-φ)*λ, (1-d)*(1-α)*(1-β), (1-d)*λ*(1-σ)],
[(1-d)*φ*(1-λ), (1-d)*β*(1-α), (1-d)*(1-σ)*(1-γ)]])
self.A_hat = self.A / (1 + self.g)
def rate_steady_state(self, tol=1e-6):
"""
Finds the steady state of the system :math:`x_{t+1} = \hat A x_{t}`
Returns
--------
xbar : steady state vector of employment and unemployment rates
"""
x = np.full(3,1/3)
error = tol + 1
while error > tol:
new_x = self.A_hat @ x
error = np.max(np.abs(new_x - x))
x = new_x
return x
def simulate_stock_path(self, X0, T):
"""
Simulates the sequence of Employment, Unemployment and Employment under TU stocks
Parameters
------------
X0 : array
Contains initial values (E0, U0, ETU0
T : int
Number of periods to simulate
Returns
---------
X : iterator
Contains sequence of employment, unemployment and employment under trade union stocks
"""
X = np.atleast_1d(X0) # Recast as array just in case
for t in range(T):
yield X
X = self.A @ X
def simulate_rate_path(self, x0, T):
"""
Simulates the sequence of employment and unemployment rates
Parameters
------------
x0 : array
Contains initial values (e0,u0)
T : int
Number of periods to simulate
Returns
---------
x : iterator
Contains sequence of employment and unemployment rates
"""
x = np.atleast_1d(x0) # Recast as array just in case
for t in range(T):
yield x
x = self.A_hat @ x
nlm = NewLakeModel()
e, f, c = np.linalg.eigvals(nlm.A_hat)
abs(e), abs(f), abs(c)
nlm = NewLakeModel()
N_0 = 200
e_0 = 0.4
eut_0 = 0.25
u_0 = 0.35
T = 50
U_0 = u_0 * N_0
E_0 = e_0 * N_0
EUT_0 = eut_0 * N_0
fig, axes = plt.subplots(4, 1, figsize=(10, 8))
X_0 = (U_0, E_0, EUT_0)
X_path = np.vstack(tuple(nlm.simulate_stock_path(X_0, T)))
axes[0].plot(X_path[:, 0], lw=2)
axes[0].set_title('Unemployment')
axes[1].plot(X_path[:, 1], lw=2)
axes[1].set_title('Employment')
axes[2].plot(X_path[:, 2], lw=2)
axes[2].set_title('Employed under Trade Union')
axes[3].plot(X_path.sum(1), lw=2)
axes[3].set_title('Labor force')
for ax in axes:
ax.grid()
plt.tight_layout()
plt.show()
nlm = NewLakeModel()
e_0 = 0.4
eut_0 = 0.25
u_0 = 0.35
T = 60 # Simulation length
xbar = nlm.rate_steady_state()
fig, axes = plt.subplots(3, 1, figsize=(10, 8))
x_0 = (u_0, e_0, eut_0)
x_path = np.vstack(tuple(nlm.simulate_rate_path(x_0, T)))
titles = ['Unemployment rate', 'Employment rate', 'Employment rate under Trade Union']
for i, title in enumerate(titles):
axes[i].plot(x_path[:, i], lw=2, alpha=0.5)
axes[i].hlines(xbar[i], 0, T, 'r', '--')
axes[i].set_title(title)
axes[i].grid()
plt.tight_layout()
plt.show()