Where did I go wrong in computing $P\left(Y-X < \frac{1}{8}\right)$ where $X \sim \mathcal{U}[0,1]$, $Y|X \sim \mathcal{U}[X,1]$? I'm interested in an adaptation of Problem 1.3 from Mathematical Statistics by S.S. Wilks (second edition, 1962). I've added the self-study tag, however this is for my amusement rather than 'homework'.

A store opens at 9 A.M. and closes at 5 P.M. A shopper taken "at random" walks into this store at time $x$ and out at time $y$ being measured in hours on the time axis with 9 A.M. as origin). Describe the sample space in terms of $(x,y)$. Describe, in terms of $x$ and $y$, the following events
(a) "The shopper is in the store less than one hour."
(b) "The shopper is in the store at time $z$."
(c) "The shopper went into the store before time $u$ and out after time $v$."

Let's focus on part (a) for this question. Rather than describe the sample spaces (in terms of set-builder notation) as the question requests, I am interested in the probability of the duration being less than 1 hour. The duration of 9 AM to 5 PM is 8 hours, so one hour is $\frac{1}{8}$th of the interval. So let us instead consider the interval $[0,1]$ and whether $y-x < \frac{1}{8}$ where $x,y \in [0,1]$ and $x < y$. Given the indicator function $\chi (x,y) = \begin{cases} 1 & y - x < \frac{1}{8} \\ 0 & \text{else} \end{cases}$ I had thought that computing $\int \int \chi (x,y) f_{X,Y|X <Y}(x,y)\ dydx$ would give me the probability of $x-y$ being less than an eighth of the interval in length.
Let the prospective density of $X$ be $f_X(x) = \frac{1}{1-0} = 1$ and and the prospective conditional density of $Y$ be $f_{Y|X}(x)=\frac{1}{1-x}$ so that $f(x,y) = f(y|x) f(x)$ is their prospective joint density function. For $f(x,y)$ to be a probability density function it must satisfy $\int_0^1 \int_0^1 f(x,y)\ dydx = 1$.
$$\begin{align}\int_0^1 \int_0^1 f(x,y)\ dydx &= \int_0^1 \int_0^1 f(y|x) f(x)\ dydx \\
 &= \int_0^1 \int_0^1 \frac{1}{1-x} \ dydx \\ &= \int_0^1 \left[ \frac{y}{1-x} \right]_{y=0}^{y=1} \ dx  \\ &= \int_0^1\frac{1}{1-x} \ dx  \\ &= -\ln |1-1| \\ &\neq 1\end{align}$$
What is clear from this is that $f(x,y)$ is not a density function inspite of the well-intentioned use of conditional probabilities. That seems to be due to the a vertical asymptote at $x=1$, which with some thought could have been noticed even before crunching any algebra.
Realizing that I might have made a mistake, possibly integration bounds or not considering a limit of the integrated area in some way, I turned to Monte Carlo methods for some intuition and sanity checks on the behaviour of this distribution. I apologize in advance for the quality of the plots. I will explain the layout of the panels and the meaning of the axes in text. The top right panel is the histogram representing $10^4$ values of $y-x$ sampled first from $x\sim \mathcal{U}[0,1]$ and then $y\sim \mathcal{U}[x,1]$. The bottom right panel represents the same, but it is the cumulative histogram. These show us intuitively (in a frequentist way) that there are more ways for shorter subintervals to be sampled than there are for longer ones to be sampled. One can compute the counting probability of $y-x < \frac{1}{8}$ from this simulation, which I did for $10^4$ iterations (thus $10^8$ resamplings from the distribution total) to get an estimate of the sampling distribution for this probability. Thus, the upper left panel is the histogram of such estimated probabilities, and the lower left panel is the cumulative histogram corresponding to the same resampled estimates. It is likely that the probability is close to 0.385, and very likely to be in $[0.37, 0.40]$.

In case of coding errors, or to improve understanding of the problem, here is the code I whipped up for the Monte Carlo method.
import numpy as np
import matplotlib.pyplot as plt


N = 10000
estimates = []
for j in range(N):
    delta = []
    for i in range(N):
        x = np.random.uniform(low=0, high=1, size=1)[0]
        y = np.random.uniform(low=x, high=1, size=1)[0]
        delta.append(y-x)
    estimates.append(np.mean(np.array(delta) < 1/8))
    print(j, estimates[-1])

fig, axes = plt.subplots(2, 2)

axes[0,0].hist(estimates, bins=min(N, 100))
axes[1,0].hist(estimates, bins=min(N, 100), cumulative=True)
axes[0,1].hist(delta, bins=min(N, 100))
axes[1,1].hist(delta, bins=min(N, 100), cumulative=True)
plt.savefig('testplot.pdf')

The Monte Carlo simulations suggest that there is some 'nice' behaviour where the discontinuity may not be important. I may have misformulated the mathematical problem and that some treatment of limits or boundaries should set me straight. Where did I go wrong?
Edit
As suggested by @whuber I drew a picture. In particular, I drew the unit square where yellow met the criterion of $x<y$ and $y-x < \frac{1}{8}$. Marginalized across all values of $(X,Y)$ we have $1 - \frac{1}{2} - \frac{(1-\frac{1}{8})^2}{2} \approx 0.12$ while conditioned on $x<y$ we have @shang Zhang's answer of $\approx 0.23$. Note that the y-axis is inverted due to how images are displayed with matplotlib.pyplot.imshow.

A remaining confusion is why both of these geometric answers disagree with the Monte Carlo simulation. Unless there is an issue in the code itself, there should be an interpretation of what that probability is and why it is different from both of the geometrical answers. Something is still missing in my understanding of the problem formulation (I suspect specifying the event carefully enough).
This confusion led me to converting the Monte Carlo simulation into the geometric picture. The same procedure as previously was followed, except for only $10^5$ and the $(x,y)$ values were plotted instead of being used to estimate the probability. The colours indicate whether the condition of $y-x < \frac{1}{8}$.

The scatter plot alone might suggest a similar answer to the geometric one given by @shang zhang, however in this Monte Carlo the KDE plots highlight an observation I made previously: not all interval lengths are equally likely. This is evidently the case for intervals that are also less than one eighth in length. Indeed, not all pairs of $(X,Y)$ meeting the conditions are equally likely.
 A: First, the joint distribution of $(X,Y)$ is indeed given by its density
$$f_X(x)\times f_{Y|X=x}(y)$$This is a general identity assuming that the distribution of the pair $(X,Y)$ is absolutely continuous with respect to the Lebesgue measure on $\mathbb R^2$ (counterexample: $Y=2X$).
Second, the error stems from missing the indicator functions in both $f_X$ and $f_{Y|X=x}$. They should be
$$f_X(x)=\mathbb I_{(0,1)}(x)\qquad f_{Y|X=x}(y)=\frac{\mathbb I_{(x,1)}(y)}{1-x}$$
under the uniform assumptions. (The phrasing of the question is too vague in this respect and another model could be proposed. For instance, Shang Zhang's answer assumes $(X,Y)$ is uniform over the upper triangle, which amounts to assume that $(X,Y)=(U_{(1)},U_{(2)})$ is the order statistic of a pair of uniforms.)

To engage with a question in the comments, the double integral
$$\int_0^1 \int_x^1 \frac{1}{1-x}~\text dy~\text dx = 1 $$ is correct
as it writes$$\int_0^1 \underbrace{\left\{\int_x^1 \frac{1}{1-x}~\text
  dy\right\}}_\text{function of $x$}~\text dx $$
meaning that the inner integral is computed for each value of $x\in(0,1)$ and then that the resulting function (of $x$) is integrated over $(0,1)$ in $x$. It is further equal to the double integral over the domain $$\{(x,y); ~0<x<y<1\}$$ as the value of an integral is independent of the order of integrands. There is thus no
ambiguity with this expression, contrary to the one in the link.

Third, for this model, the value of $$\mathbb P(Y-X<1/8)$$can be found by

*

*integrating out the indicator function $\mathbb I_{x,\min\{1,x+1/8\}}(y)$ with respect to the joint density, i.e.,
\begin{align}\int_0^1\int_x^{\min\{1,x+1/8\}} \frac{\mathbb I_{(x,1)}(y)}{1-x}\,\text d y~\text dx&=\int_0^1 \frac{\min\{1,x+1/8\}-x}{1-x}~\text dx\\
&=\frac{1}{8}-\frac{1}{8}\log\frac{1}{8}
\end{align}

*deriving the distribution of $Z=Y-X$
$$Z\sim -\log(z)\mathbb I_{0<z<1}$$
and computing$$\mathbb P(Z<1/8)=-\int_0^{\frac{1}{8}}\log z\text dz=\frac{1}{8}-\frac{1}{8}\log\frac{1}{8}$$
A: Split in two cases:
$
P(Y-X \leq z) = \int_0^1 P(Y \leq z + x \mid X = x)f_X(x)=\int_0^1 P(A(x))\cdot 1\\=\int_0^{1-z} P(A(x))+ \int_{1-z}^1 P(A(x)) = I_1 + I_2.
$
$Y$ is smaller than 1 so $z+x>1$ implies that $P(A(x)) = 1$
and $I_2 = \int_{1-z}^1 1 \,dx = z$.
$P(Y \leq z + x \mid X = x) = \frac{z}{1-x}\int_0^{1-z}dy= \frac{z}{1-x}$ so
$I_1 = z\int_0^{1-z}\frac{1}{1-x} = -z\log(z)$.
Therefore
$P(Y-X \leq z) = I_1 + I_2 = z(1-\log(z))$.
$\frac{1}{8} + \frac{1}{8}\log{8}\approx 0.384930$
A: I apologize for not reading through your entire post (as it is pretty long). But I think I can give you the answer to the problem: "Draw $X$ uniformly from [0, 1], and then, conditional on that $Y$ uniformly from $[X, 1]$, what is the probability that they are within 1/8?"
Answer: Imagine you randomly pick a point from the upper triangle of the 1-by-1 square. (This is the data generating process for your $(X, Y)$. The area of the triangle is 1/2. If this point is NOT within 1/8 of the diagonal line, then it is in the upper left corner triangle with area 1/2 * 7/8 * 7/8 = 49/128.
So the answer = 1 - (49/128) / (1/2) = 15/64.
