Assume that the waiting time for each event follows the exponential distribution with the parameter $\lambda$. Let $\tau_{i,j}$ denote the waiting time between $i$th event and $j$th event. For a fixed duration $T$, what is the probability $P(\tau_{1,2}+\tau_{2,3}<T \text{ and }\tau_{1,2}+\tau_{2,3}+\tau_{3,4}>T)$? That is, Events 2 and 3 happened within time T after the occurrence of Event 1, but Event 4 has not.
I am stuck in finding the joint pdf of $\tau_{1,2}+\tau_{2,3}<T$, and $\tau_{1,2}+\tau_{2,3}+\tau_{3,4}>T$.
Let $X=\tau_{1,2}+\tau_{2,3}$ and $Y=\tau_{3,4}$, then its equivalent to find $P(X<T \text{ and } Y>T-X)$, here $X$ is a random variable follows the gamma distribution.
If I have the joint pdf, then I can compute $P(X<T \text{ and } Y>T-X)$ as follows: $\int_{0}^{T}\int_{T-x}^{\infty} f_{X,Y}(x,y) \,dy\,dx$.
How to find the $f_{X,Y}(x,y)$ in this case?
Assuming independence of waiting times, the joint density for three random variables $(\tau_{1,2},\tau_{2,3},\tau_{3,4})$ when they are positive is the product of their individual densities $$f(\tau_{1,2},\tau_{2,3},\tau_{3,4})=\lambda^3 e^{-\lambda(\tau_{1,2}+\tau_{2,3}+\tau_{3,4})}.$$
The tricky part is setting the limits of integration. I think you are looking for
$$\begin{align}&P(\tau_{1,2}+\tau_{2,3}<T \text{ and }\tau_{1,2}+\tau_{2,3}+\tau_{3,4}>T) \\=&\int_{\tau_{1,2}=0}^T\int_{\tau_{2,3}=0}^{T-\tau_{1,2}}\int_{\tau_{3,4}=T-\tau_{1,2}-\tau_{2,3}}^{\infty} f(\tau_{1,2},\tau_{2,3},\tau_{3,4})\, d\tau_{3,4}\,d\tau_{2,3}\,d\tau_{1,2}\\ =& \frac12(\lambda T)^2 e^{-\lambda T}\end{align}$$
This is a probability, but happens to be proportional to the expression for the density for a gamma or Erlang distribution of shape $3$ (and identical when $\lambda=1$).
Alternatively you could see this as the difference between two CDFs of gamma or Erlang distributions with different results: you want $P(\tau_{1,2}+\tau_{2,3}\le T) - P(\tau_{1,2}+\tau_{2,3}+\tau_{3,4}\le T)$, with subtraction possible since the latter can only happen when the former does.
As you would hope, these provide the same result. Plotting the two curves with R gives a perfect match.