Solving Kaplan and Menzio: Shopping Time

Question

Kaplan and Menzio's shopping time model is a search&matching unemployment model where we, for a steady state equilibrium, need to determine to variables:

• $J$: Value of a worker
• $u$: unemployment rate

Their steady state values ($\dot J, \dot u = 0$) are given by the following system of equations:

$$ J = \frac{(1-\gamma)(S(u) + y_e - y_u)}{\rho + \delta}\\ u = \frac{\delta}{\delta + \lambda(J)} $$

Furthermore, "profits from the good with price dispersion" and wages, respectively, are given by (equations (6), (7) on page 12)

$$ S(u) = A(u) + B(u)*w(u)\\ w(u) = y_u + \gamma*(S(u) + y_e - y_u) $$

Where $A(u), B(u)$ are some placeholders for a bunch of constants and parameters, that may depend on $u$. I plug out the wages and solve as

$$ \tilde S(u) = \frac{A(u) + B(u)(y_u + \gamma(y_e - y_u)}{1 - B(u)\gamma} $$

Then, using the authors' calibration and this version of $\tilde S(u)$, I can solve for both $J$ and $u$. Here it is:

Given their calibration, they get two interesting steady state equilibria on the $(u, J)$ space: $(0.053, 11.3)$ and $(0.081, 9.5)$. My $u(J)$ graph seems to be correct as it goes through both these loci. The $J(u)$ line, however, seems to be "too small" given $u$.

The only object somewhat difficult to compute is $S(u)$. I have verified that my step of replacing actually works: I computed $w(u)$ using $\tilde S(u)$ and verified that $S(u, w(u)) == \tilde S(u)$ (up to a numerical approximation).

I am not asking for someone to find my mistake, but rather where to look for it: I have checked anything I could come up with and am somewhat stuck. For the interested, here is my Python code:

# solves simple shopping time model
# only on the steady state equilibria

import numpy as np
import matplotlib.pyplot as plt

# all parameters for monthly
class Parameters(object):
    def __init__(self):

        self.psi_u = 0.27 # search twice, unemployed
        self.psi_e = 0.02 # search twice, employed < psi_u
        self.gamma = 0.74 # bargaining power
        self.alpha = 1 # BJ exponent in utility

        self.r = 15.7 # reservation price, home transformation
        self.c = 1 # production technology firm

        self.y_u = 4.91 # efficiency home
        self.y_e = 1 # efficiency firm

        self.rho = 0.003 # depreciation rate
        self.delta = 0.024 # separation rate

        self.phi = 1.24 # matching function efficiency
        self.k = 8.02 # search costs of labor for firm


def b(u, Param):
    return 1 + Param.psi_e + u*(Param.psi_u - Param.psi_e)
def s(u, Param):
    return 1 - u
def sigma(u, Param):
    return s(u, Param)/b(u, Param)
# seller meets buyer
def muOfU(u, Param):
    return np.minimum(1/sigma(u, Param), 1)
# buyer meets seller
def nuOfU(u, Param):
    return np.minimum(sigma(u, Param), 1)


def lambdaOfJ(J, Param):
    result = ( ((J/Param.k)**(Param.phi) - 1 )**(-1) + 1 )**(-float(1)/Param.phi)
    result[np.isnan(result)] = 0
    return result

# not needed atm, replaced within lambdaOfJ (instead of lambdaOfTheta)
def thetaOfJ(J, Param):
    #eta(theta)J = k
    #etaInverse = (eta^(-phi) - 1)^(1/phi)
    x = Param.k/J
    return (x**(-Param.phi) - 1)**(1/Param.phi)



def STilde(u, Param):


    A1 = muOfU(u, Param)*u*(1+Param.psi_u)/b(u, Param)
    A2 = 1 - 2*Param.psi_u*nuOfU(u, Param)/(1+Param.psi_u)
    A3 = Param.alpha*Param.y_u*(Param.r-Param.c)/Param.r
    A = A1*A2*A3

    B1 = muOfU(u, Param)*(1-u)*(1+Param.psi_e)/b(u, Param)
    B2 = 1 - 2*Param.psi_e*nuOfU(u, Param)/(1+Param.psi_e)
    B3 = Param.alpha*(Param.r-Param.c)/Param.r

    B = B1*B2*B3

    result = (A + B*(Param.y_u + Param.gamma*(Param.y_e - Param.y_u) ))/(1 - B*Param.gamma)
    return result

def S(u, wages, Param):
    A1 = muOfU(u, Param)*u*(1+Param.psi_u)/b(u, Param)
    A2 = 1 - 2*Param.psi_u*nuOfU(u, Param)/(1+Param.psi_u)
    A3 = Param.alpha*Param.y_u*(Param.r-Param.c)/Param.r
    A = A1*A2*A3

    B1 = muOfU(u, Param)*(1-u)*(1+Param.psi_e)/b(u, Param)
    B2 = 1 - 2*Param.psi_e*nuOfU(u, Param)/(1+Param.psi_e)
    B3 = Param.alpha*(Param.r-Param.c)/Param.r

    B = B1*B2*B3

    result = A + B*wages
    return result



def w(u, Param):
    return Param.y_u + Param.gamma*(STilde(u, Param) + Param.y_e - Param.y_u)

#result2 = A + B*w(u, Param)

def JSteadyState(u, Param):

    A = (1-Param.gamma)*(STilde(u, Param) + Param.y_e - Param.y_u )
    B = Param.rho + Param.delta
    return A/B

def uSteadyState(J, Param):

    B = Param.delta + lambdaOfJ(J, Param)
    return Param.delta/B


Param = Parameters()


# compute J(u) and u(J)
uGrid = np.linspace(0.001, 0.2, 100)
J = JSteadyState(uGrid, Param)

JGrid = np.linspace(5, 20, 1000)
u = uSteadyState(JGrid, Param)


# verify computation of S(u)
# wages given STilde:
wages = w(uGrid, Param)
# S given u, STilde
SOriginal = S(uGrid, wages, Param)
SWithTilde = STilde(uGrid, Param)

print max(abs(SWithTilde - SOriginal)) # 8.881784197e-16

# plot J(u) and u(J)
plt.plot(u, JGrid, label='u(J)')
plt.plot(uGrid, J, label='J(u)')
plt.legend()
plt.xlabel('u')
plt.ylabel('J')
plt.show()

Answer 1

UPDATE
After e-mail communication with one of the authors G.W.Kaplan, I recalibrated the value of the vacancy-posting cost parameter $k$ in order to obtain a cross of the two nullclines for $u=0.05$. This is achieved for $k=7.41$ (rounded). Moreover, with this value of $k$, I get a second (but not a third) steady state. A close up diagram :

This still is not what Figures 1 and 2 of the paper show. The shape of the $J$-nullcline remains the issue. It should have an initial concave part and then a convex part in order to give us three steady states. But I won't dig it any further than that.

---

The $u$-nullcline

On p.8 the authors define $\lambda(\theta) \equiv M(1,\theta)$ (I believe they should have written for clarity, $\lambda(\theta) \equiv M/u = M(1,\theta)$) . On p. 25 they specify the matching function as $M(u,v) = uv(u^\phi + v^\phi)^{-1/\phi}$. Therefore we have

$$\lambda(\theta) = \theta(1 + \theta^\phi)^{-1/\phi}$$

This means that whenever $\theta = 0 \implies \lambda(\theta) =0$.

On p.14, eq. $(8)$ the give $k> J \implies \theta = 0$. In the calibration stage, $k=8.02$. The $u$-nullcline is defined as

$$u_{SS} = \frac {\delta}{\delta + \lambda[\theta(J))]}$$

All these together imply that $$J < 8.02 \implies u_{SS} = 1$$

Using the other relations specified, after a little algebra, we get $$J>k=8.02 \implies \lambda[\theta(J)] = \left(1-(k/J)^{\phi}\right)^{1/\phi}$$

which is indeed your code, a bit less convoluted. And indeed the $u$-nullcline never crosses with the $J$-nullcline. For values of $u$ near unity see below. For low values of $u$ we obtain

\begin{array}{| r | r |} \hline \hline \text {u} & \text {Jss} & \text{uSS(in J units)} \\ \hline 0.03 & 12.25 & 23.06 \\ 0.04 & 11.53 & 14.13 \\ 0.05 & 10.86 & 11.75 \\ 0.06 & 10.23 & 10.66 \\ 0.07 & 9.64 & 10.03\\ 0.08 & 9.09 & 9.62\\ 0.09 & 8.57 & 9.34\\ 0.10 & 8.09 & 9.14\\ \hline \end{array}

The $u$-nullcline stays always above the $J$-nullcline. I also "pushed" the $J$-nullcline by the factor $(1+\rho)$, i.e. the discrete version of the equation. But $\rho=0.003$, too small a mark-up, and so again there was no crossing. Therefore it appears that what went into the simulations of the authors is at some point(s) different from what the equations and the calibrated parameters (Table 2) appearing in the paper give us. Or, we are missing something.

The $J$-nullcline
I checked the OP's code line by line, including the values of the calibrated parameters. I did not find any discrpeancy with what the paper gives. I then copy-pasted the code for the $J$-nullcline into Gretl, tweaking it only to match the local language. A note: the way various magnitudes are defined in the paper, we have

$$\psi_u > \psi_e \implies b >1 \;\forall u \implies \sigma = s/b < 1, \;\forall u$$

$$\implies \text {muOfU}\equiv \mu(\sigma(u)) = \min(1/\sigma, 1) = 1,\; \forall u$$

$$\implies \text {nuOfU}\equiv \nu(\sigma(u)) = \min(\sigma, 1) = \sigma, \; \forall u$$

**GRETL script**  
solves simple shopping time model
 only on the steady state equilibria
     nulldata 100  
      genr series u = index/100   
        genr scalar  psi_u = 0.27 # search twice, unemployed
        genr scalar  psi_e = 0.02 # search twice, employed < psi_u
        genr scalar  gamma = 0.74 # bargaining power
        genr scalar  alpha = 1 # BJ exponent in utility

        genr scalar  r = 15.7 # reservation price, home transformation
        genr scalar  c = 1 # production technology firm

        genr scalar  y_u = 4.91 # efficiency home
        genr scalar  y_e = 1 # efficiency firm

        genr scalar  rho = 0.003 # depreciation rate
        genr scalar  delta = 0.024 # separation rate

        genr scalar  phi = 1.24 # matching function efficiency
        genr scalar  k = 8.02 # search costs of labor for firm

        genr series b =1 + psi_e + u*(psi_u - psi_e)
        genr series s = 1-u 
        genr series sigma = s/b
        genr series sigmareci = 1/sigma
        genr list  mlist = sigmareci const
        genr list  nlist  = sigma const
        genr series muOfU = min(mlist)  # seller meets buyer
        genr series nuOfU = min(nlist)  # buyer meets seller


# def STilde:

    genr series A1 = muOfU*u*(1+psi_u)/b
    genr series A2 = 1 - 2*psi_u*nuOfU/(1+psi_u)
    genr scalar A3 = alpha*y_u*(r-c)/r
    genr series A = A1*A2*A3

    genr series B1 = muOfU*(1-u)*(1+psi_e)/b
    genr series B2 = 1 - 2*psi_e*nuOfU/(1+psi_e)
    genr series B3 = alpha*(r-c)/r

    genr series  B = B1*B2*B3

    genr series STilde = (A + B*(y_u + gamma*(y_e - y_u) ))/(1 - B*gamma)

# def JSS

genr series JSS = (1-gamma)*(STilde + y_e - y_u ) / (rho+delta)

I obtained the values

and the graph

From the above we can see that the $J$-nullcline, although it starts increasing, it never approaches again the value $8$, which is $\approx$ the low asymptote for the $u$-nullcline.

Answer 2

It is possible that this issue is caused by the limits of the precision of the calculation operations. The Python package bigfloat may be of interest in ruling this out.