What is the statistical steady state of this poisson process? 
Suppose there is an island with some amount of people. The arrival of new people on the island follows a Poisson process with rate parameter $\lambda$. Suppose further that the time until a person departs the island is exponential with rate constant $\mu$. Assume all arrival and departure times are independent. An event is the arrival of a new person or departure of an existing person.
a) Suppose there are $n$ people on the island at time $t$. What is the probability distribution of time until the next event?
b) What is the probability that the next event is a departure?
c) Assume a statistical steady state exists. What is the expected number of people  in this steady state?

Attempt
Assuming my attempt on a) and b) is correct, I would like to understand better how to approach c). 
For a), let $T$ be the time until the next event. Then
$$ P(T<s) = 1 - e^{-s(\lambda + \mu n)}$$
so the density is given by
$$ (\lambda + \mu n)e^{-s(\lambda + \mu n)}.$$
Hence $T$ is exponential with rate constant $(\lambda + \mu n)$.
For b), I discretized time and examined
$$ P( T = k \, \Delta s \cap N_{t+k\,\Delta s} = N_t ) $$
where $N_{t+k\,\Delta s} = N_t$ is the number of new people in time $k \,\Delta s$. This leads to
$$ P( \text{Next event is a departure} ) = 
\frac{\lambda + \mu n }{ 2\lambda + \mu n}.$$
Which seems to make sense, since for large $n$, the probability that the next event is a departure is large.
Now, for part c), I'm not sure how to marry the two concepts above to arrive at a formula for the statistical steady state. I expect I need to arrive a formula for a time $t$ and take $t \to \infty$. If my solutions for a) and b) are right, is there a logical way to use them to produce the answer for c)?
 A: Let $\{N(t):t\geqslant0\}$ be the number of people on the island at time $t$, and $\{X_m:m=1,2,\ldots\}$ be time of the $m^{\mathrm{th}}$ event. 
Since the arrival and departure processes are Markovian, $N(t)$ is a continuous-time Markov chain with embedded Markov chain $X_m$, having transition probabilities
$$
\mathbb P(X_{m+1}=j\mid X_m=i) = \begin{cases}
\frac\lambda{\lambda+i\mu},& j=i+1\\
\frac{i\mu}{\lambda+i\mu},& j=i-1.
\end{cases}
$$
So given $N(t)=n$, it follows that
$$\mathbb P(T>s) = e^{-(\lambda+n\mu)s},\ s>0. $$
Let $q_{ij}$ be the transition rate from state $i$ to state $j$. If a stationary distribution $\pi$ exists, then $\pi$ must satisfy the global balance equations $$\sum_{j\in S\setminus\{i\}}\pi_i q_{ij} = \sum_{j\in S\setminus\{i\}}\pi_j q_{ji} $$
for each $j\in S$ (here the state space is $S=\{0,1,2,\ldots\}$). Since
$$
q_{ij} = \begin{cases}
\lambda,& i=j+1\\
(i+1)\mu,& i=j-1,
\end{cases}
$$
the balance equations are given by
$$\lambda\pi_n = (n+1)\mu \pi_{n+1},\ n\geqslant 0, $$
which yields the recurrence
$$\pi_n =\left(\frac\lambda\mu\right)^n\left(\frac1{n!}\right)\pi_0,\ 
n\geqslant0. $$
Set $\rho = \frac\lambda\mu$. From $\sum_{n=0}^\infty \pi_n=1$ it follows that 
$$\pi_0\sum_{n=0}^\infty\frac{\rho^n}{n!}=1, $$
and hence 
$$\pi_n = \frac{e^{-\rho}\rho^n}{n!}. $$
By inspection, $\pi$ is a Poisson distribution with parameter $\rho$, and so $$\lim_{t\to\infty}\mathbb E[N(t)] = \rho. $$
Alternatively, we can compute the stationary distribution directly from the transition rates. For each $i,j\geqslant 0$, define the generator (or transition rate matrix) by
$$
G_{ij} = \begin{cases}
q_{ij},& i\ne j\\
-\sum_{j\ne i}q_{ij},& i=j
\end{cases}
$$ 
Let $P(t)$ be the matrix with entries $\mathbb P(N(t)=j\mid N(0)=0)$.Then $P$ satisfies the differential equations
\begin{align}
[P'(t)]_{ij}&=[P(t)G]_{ij}, P(0) = \mathsf I,\quad t>0\\
[P'(t)]_{ij}&=[GP(t)]_{ij}, P(0) = \mathsf I,\quad t>0
\end{align}
i.e. the Kolmogorov forward and backward equations. The solution is given by
$$P(t) = e^{tG} = \sum_{n=0}^\infty \frac{(tG)^n}{n!}, $$
and so for any $j\in S$, $$\lim_{t\to\infty}P(t)_{ij} = \pi_j $$ is the stationary distribution of $N(t)$.
Now, to actually compute $\pi$, we use the following result:
$$\pi = \pi P \iff \pi G=0. $$
(This is a good exercise to prove.). Here the generator matrix is given by
$$G_{ij} = \begin{cases}
\lambda,& i=j-1
-(\lambda+i\mu),& i=j 
(i+1)\mu,& i=j+1,\end{cases}
$$
so we have
\begin{align}
-\lambda\pi_0 + \mu\pi_1=0&\implies \pi_1 = \rho\pi_0,\\
\lambda\pi_0 - (\lambda+\mu)\pi_1 + 2\mu\pi_2&\implies \pi_2 = \frac{\rho^2}{2!}\pi_0,
\end{align}
and for $n\geqslant 2$, induction yields
$$\pi_n = \frac{\rho^n}{n!}\pi_0. $$
As before,
$$\pi_0\sum_{n=0}^\infty\frac{\rho^n}{n!}=1\implies \pi_0 = e^{-\rho}  $$
and hence
$$\pi_n = \frac{e^{-\rho}\rho^n}{n!}. $$
