Consider three independent Poisson processes $N_1(t), N_2(t),N_3(t)$ with rates $\lambda_1,\lambda_2,\lambda_3$, respectively. What is the probability of reaching $(i,j,k)$?
I know that for the two-dimensional process, we have the probability of reaching $(i,j)$ as $$\frac{(i+j)!\lambda_1^i\lambda_2^j}{i!j!(\lambda_1+\lambda_2)^{i+j}}$$ Could the problem be easily extended to the three-dimensional Poisson process as $$\frac{(i+j+k)!\lambda_1^i\lambda_2^j\lambda_3^k}{i!j!k!(\lambda_1+\lambda_2+\lambda_3)^{i+j+k}}$$
I have reason to believe so since by considering $\lambda_1=\lambda_2=\lambda_3=1$; $i=j=k=n$, we have the probability of reaching $(n,n,n)$ as $$\frac{3n!\left(\frac{1}{3}\right)^{3n}}{n!n!n!}$$
$$\int_0^t P(N_1(t)=i,N_2(t)=j,N_3(t)=k) dt$$ is the average amount of time the process spends in $(i,j,k)$. This is equal to:
$$(*)\quad P(\mbox{ever reaching }(i,j,k)) \times \frac{1}{\lambda_1+\lambda_2+\lambda_3},$$
because time process spends at any state it arrives to is $\mbox{Exp}(\lambda_1+\lambda_2+\lambda_3)$.
This is also equal to
\begin{align*} &\int_0^\infty P(N_1(t) =i) P(N_2(t)=j) P(N_3(t)=k) dt\\ & =\int_0^\infty e^{-(\lambda_1+\lambda_2+ \lambda_3)t} t^{i+j+k} \frac{ \lambda_1^i\lambda_2^j \lambda_3^k}{i! j! k!}d t\\ &= \frac{ \lambda_1^i\lambda_2^j \lambda_3^k}{i! j! k!}(-1)^{i+j+k} \frac{\partial ^{i+j+k}}{\partial \lambda^{i+j+k}}\int_0^\infty e^{-(\lambda_1+\lambda_2+\lambda_3)t}dt \\ (**) \quad& = \frac{ \lambda_1^i\lambda_2^j \lambda_3^k}{i! j! k!}(i+j+k)!(\lambda_1+\lambda_2+\lambda_3)^{-(i+j+k+1)}. \end{align*}
Equating $(*)$ to $(**)$ we obtain
$$P(\mbox{ever reaching }(i,j,k))= \frac{(i+j+k)!}{i! j! k!}\frac{\lambda_1^i\lambda_2^j\lambda_3^k}{(\lambda_1+\lambda_2+\lambda_3)^{i+j+k}}.$$