<p>It is possible that this issue is caused by the limits of the precision of the calculation operations. The Python package <a href="https://pythonhosted.org/bigfloat/" rel="nofollow">bigfloat</a> may be of interest in ruling this out. </p>

CC BY-SA 3.0
2015-03-22T22:04:04.787Z
2015-03-22T22:04:04.787Z
BKay
<p><strong>UPDATE</strong><br>
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 :</p>

<p><img src="https://i.sstatic.net/kWwjw.png" alt="enter image description here"></p>

<p>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 <em>three</em> steady states. But I won't dig it any further than that.</p>

<hr>

<p><strong>The $u$-nullcline</strong>  </p>

<p>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</p>

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

<p>This means that whenever $\theta = 0 \implies \lambda(\theta) =0$.</p>

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

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

<p>All these together imply that
$$J &lt; 8.02 \implies u_{SS} = 1$$</p>

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

<p>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</p>

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

<p>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. 
<strong>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.</strong> Or, we are missing something.</p>

<p><strong>The $J$-nullcline</strong><br>
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.
<strong>A note:</strong> the way various magnitudes are defined in the paper, we have</p>

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

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

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

<pre><code>**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 &lt; 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)
</code></pre>

<p>I obtained the values</p>

<p><img src="https://i.sstatic.net/QedFo.png" alt="enter image description here"></p>

<p>and the graph</p>

<p><img src="https://i.sstatic.net/JIlSb.png" alt="enter image description here"></p>

<p>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. </p>

CC BY-SA 3.0
2015-03-22T21:12:31.000Z
Ad 1): Actually, the authors say at some point *In the remainder of the paper [...] it is more convenient to work in continuous time.*. Ad 2): Yes, I expect that to be the problem: But what can possibly go wrong? I show that $S(u, w(\tilde S, u)) == \tilde S(u)$, which to me indicates that $S(u, w)$ is computed correctly - and besides the $S$ term, there's not really much more to that equation.
FooBar
CC BY-SA 3.0
2015-03-22T23:56:54.550Z
About 1), I did read that quote, but sometimes such comments do not extend to actual simulations -they may meant only the theoretical explorations. As for 2) I meant more of a shallow mistake like a sign, or things like that in the actual script you wrote -not some structural mistake.
Alecos Papadopoulos
CC BY-SA 3.0
2015-03-28T17:13:51.490Z
Yes, your assertion is correct. You get exactly the same values as I do - but these are too small: We don't hit either `(0.081, 9.5)`, nor `(0.05, 11.3)` - around those two points, we instead get `J` values of `9.1` and `10.8`
FooBar
CC BY-SA 3.0
2015-03-28T18:34:40.870Z
@FooBar Then i would like also to double-check the $u$-nullcline, just in case. Can you tell me where in the paper do the authors provide explicit functional forms for the $\lambda()$ and the $\theta()$ functions? I find general definitions, not implementable  functional forms.
Alecos Papadopoulos
CC BY-SA 3.0
2015-03-28T19:03:39.750Z
You find $\lambda$ and $\theta$ as a function of $M$ on page 7. $M$ is then defined to be $M(u,v) = uv(u^\phi + v^\phi)^{-1/\phi}$ on page 24.
FooBar
CC BY-SA 3.0
2015-03-30T01:35:46.920Z
I was adding a huge bounty here, hoping that I could boost some youngsters account. Of all people, I guess you need the points the least - but you deserve them ;)
FooBar
CC BY-SA 3.0
2015-03-30T01:46:45.807Z
Thanks FooBar. You are right about the reputation points, but, indeed, I did some work here. Let's hope that it may function as some sort of example on how to deal with such questions and issues that arise in actual research. And the mystery hasn't been solved yet.
Alecos Papadopoulos
CC BY-SA 3.0
2015-03-22T18:05:21.467Z
2015-03-30T00:06:22.430Z
2015-03-30T00:06:22.430Z
Alecos Papadopoulos
<p>Kaplan and Menzio's <a href="http://se.shufe.edu.cn/upload/htmleditor/File/140522081313.pdf" rel="nofollow noreferrer">shopping time model</a> is a search&amp;matching unemployment model where we, for a steady state equilibrium, need to determine to variables:</p>

<ul>
<li>$J$: Value of a worker </li>
<li>$u$: unemployment rate</li>
</ul>

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

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

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

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

<p>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 </p>

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

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

<p><img src="https://i.sstatic.net/GgtM7.png" alt="My computation"></p>

<p>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$. </p>

<p>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). </p>

<p>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:</p>

<pre><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 &lt; 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()
</code></pre>

CC BY-SA 3.0
2015-03-22T17:28:22.580Z
The link you provide downloads a .pdf which may be problematic to view as regards mathematical notation. Google Scholar gives also the following freely downloadable link which shows no such problems, http://se.shufe.edu.cn/upload/htmleditor/File/140522081313.pdf
Alecos Papadopoulos
<macroeconomics><academic-graduate><general-equilibrium><computation><unemployment>
Solving Kaplan and Menzio: Shopping Time
FooBar
# Solving Kaplan and Menzio: Shopping Time

(Asked by: username_0 on 2015-03-20)

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()

**Tags:** macroeconomics, academic-graduate, general-equilibrium, computation, unemployment

> **Comments:**
> **username_1**: The link you provide downloads a .pdf which may be problematic to view as regards mathematical notation. Google Scholar gives also the following freely downloadable link which shows no such problems, http://se.shufe.edu.cn/upload/htmleditor/File/140522081313.pdf  

---

## 2 answers

### Answer by username_1 on 2015-03-22. Score: 6 | Accepted answer

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.

> **Comments:**
> **username_0**: Ad 1): Actually, the authors say at some point *In the remainder of the paper [...] it is more convenient to work in continuous time.*. Ad 2): Yes, I expect that to be the problem: But what can possibly go wrong? I show that $S(u, w(\tilde S, u)) == \tilde S(u)$, which to me indicates that $S(u, w)$ is computed correctly - and besides the $S$ term, there's not really much more to that equation.  
> **username_1**: About 1), I did read that quote, but sometimes such comments do not extend to actual simulations -they may meant only the theoretical explorations. As for 2) I meant more of a shallow mistake like a sign, or things like that in the actual script you wrote -not some structural mistake.  
> **username_0**: Yes, your assertion is correct. You get exactly the same values as I do - but these are too small: We don't hit either `(0.081, 9.5)`, nor `(0.05, 11.3)` - around those two points, we instead get `J` values of `9.1` and `10.8`  
> **username_1**: @username_0 Then i would like also to double-check the $u$-nullcline, just in case. Can you tell me where in the paper do the authors provide explicit functional forms for the $\lambda()$ and the $\theta()$ functions? I find general definitions, not implementable  functional forms.  
> **username_0**: You find $\lambda$ and $\theta$ as a function of $M$ on page 7. $M$ is then defined to be $M(u,v) = uv(u^\phi + v^\phi)^{-1/\phi}$ on page 24.  
> **username_0**: I was adding a huge bounty here, hoping that I could boost some youngsters account. Of all people, I guess you need the points the least - but you deserve them ;)  
> **username_1**: Thanks username_0. You are right about the reputation points, but, indeed, I did some work here. Let's hope that it may function as some sort of example on how to deal with such questions and issues that arise in actual research. And the mystery hasn't been solved yet.  

---

### Answer by username_2 on 2015-03-22. Score: 2

It is possible that this issue is caused by the limits of the precision of the calculation operations. The Python package [bigfloat](https://pythonhosted.org/bigfloat/) may be of interest in ruling this out.


---

