This is not a full answer. I found a different perspective to view the problem in, which turns the question from a complicated integral to a complicated discrete summation.
Consider a Poisson process, with unit density, on the positive real line. This means that for any $0<a<b$, the number of arrivals in $[a,b]$ has the distribution $\text{Poisson}(b-a)$. If you number the arrivals $0<A_1<A_2<\dots$, with the convention $A_0=0$, then it is well known that the inter-arrival times $A_k-A_{k-1}$ all have an exponential distribution with rate $1$. So, from now on, I will identify $$X_k=A_{k}-A_{k-1}.$$
Now, how do your events related to the Poisson process?
You are requiring $(X_1-1)\ge 0$, which implies that $A_1=X_1\ge 1$. That is, the first arrival happens after time $1$, so there are no arrivals in the interval $(0,1)$.
Next, you require $(X_1-1)+(X_2-1)\ge 0$, which implies $A_2=X_1+X_2\ge 2$. That is, the second arrival happens after time $2$, there is at most one arrival in the interval $(0,2)$.
Continuing in this fashion, for each $k\in \{1,\dots,n\}$, the event $(X_1-1)+\dots+(X_k-1)\ge 0$ is equivalent to the Poisson process having at most $k-1$ arrivals in the interval $(0,k)$. To simplify this further, for each $k\in \{1,\dots,n\}$, define
$Z_k$ to be the number of arrivals in $(k-1,k)$. Then, let
$$
E_k=\{Z_1+Z_2+\dots+Z_k\le k-1\}.
$$
We have shown that
$$
p_n=\mathbb P(E_1\cap E_2\cap \dots \cap E_n).
$$
Let us see how this plays out when $n=3$. According to above,
$$
P_3=\mathbb P(Z_1\le 0\;\cap\; Z_1+Z_2\le 1\;\cap\;Z_1+Z_2+Z_3\le 2)
$$
For these three inequalities to be satisfied, the vector $(Z_1,Z_2,Z_3)$ must be one of five possibilities: $(0,0,0),(0,0,1),(0,0,2),(0,1,0),$ or $(0,1,1)$. We can easily find the probability of each of these vectors, since $Z_1,Z_2,Z_3$ are independent, and each has a $\text{Poisson}(1)$ distribution. The result is
$$
{e^{-3}}\left(\frac1{0!\cdot0!\cdot0!}+\frac1{0!\cdot0!\cdot1!}+\frac{1}{0!\cdot0!\cdot2!}+\frac{1}{0!\cdot1!\cdot0!}+\frac{1}{0!\cdot1!\cdot1!}\right)=\frac92e^{-3}
$$
For general $n$, you need to enumerate all sequences $(z_1,\dots,z_n)$ with the property that the $z_1+\dots+z_k\le k-1$ for each $k\in \{1,\dots,n\}$, and then compute the sum of
$$
e^{-n}\frac1{z_1!\cdot z_2!\cdots z_n!}\tag1
$$
over all such valid sequences.
I find it intriguing that for each $n$, the number of valid sequences is always equal to the $n^\text{th}$ Catalan number! Specifically, we can give a bijection between valid sequences $(z_1,\dots,z_n)$, and lattice paths from $(0,0)$ to $(n,n)$ which never go above the line $y=x$. I give more details in the answer to this question.
To see the combinatorial nature of this problem, note that proving the sum $(1)$ over all valid sequences is equal to $e^{-n}n^{n-1}/(n-1)!$, is equivalent to proving
$$
{n^n}=\sum_{\substack{{z_1,\dots,z_n\ge0}\\ {z_1+\dots+z_k\le k-1\;\forall k}}}\frac{n!}{z_1!\cdots z_n!}
\tag2$$
Note that the left hand side enumerates a nice combinatorial object, namely the set of functions from an $n$-element set to itself. Furthermore, the right hand side is a sum of integers, and it also has a nice combinatorial interpretation. Therefore, I suspect that you can give a bijective proof of $(2)$. However, I cannot imagine what the bijection is.
