Finding covariance for an M/G/$\infty$ queuing system I am studying for an exam, and one of the recommended example problems is as follows:
Consider an M/G/$\infty$ queuing system. The arrivals form a Poisson Process $\mathbf{N} = \{N(t); t \ge 0\}$ having rate $\lambda >0$. The service times $S_1, S_2, \ldots$ form a sequence of iid random variables having distribution function $F$ and are independent of the arrival process. Assume that both the expected service time and the variance of a service time are strictly positive and finite. For $t \ge 0$, let $X(t)$ denote the number of servers busy at time $t$. Assume $X(0)=0$, find, for $0 < s<t$, the covariance of $X(s)$ and $X(t)$. (Note that some of the customers in the system at time $s$ will still be there at time $t$). 
My attempt
Let the entering time for a customer be $r$. This can be consider a Poisson thinning problem where a Type I customer is a customer that has completed service by time $t$, and a Type II customer is one that has not completed service by time $t$. Let $D$ be the random variable that represents the time it takes for a customer to complete service. Then one theorem we have is that 
$$\mathbb{E}[X(t)] = \lambda \int_0^t \mathbb{P}(D > t-r)dr = \lambda \int_0^t (1- F(t-r))dr$$
Similarly, 
$$\mathbb{E}[X(s)] = \lambda \int_0^s \mathbb{P}(D > s-r)dr = \lambda \int_0^s (1- F(s-r))dr$$
Then we have that 
\begin{align*}
Cov(X(s), X(t)) &= \mathbb{E}[X(s)X(t)] - \mathbb{E}[X(s)]\mathbb{E}[X(t)] \\
&= \mathbb{E} \Big[X(s)\big(X(t)-X(s)+X(s)\big)\Big] - \lambda^2 \int_0^s (1-F(s-r))dr \int_0^t (1-F(t-r))dr\\
&= \mathbb{E} [X(s)^2] + \mathbb{E}[X(s)(X(t)-X(s))] - \lambda^2 \int_0^s (1-F(s-r))dr \int_0^t (1-F(t-r))dr\\
&= Var(X(s)) + \mathbb{E}[X(s)]^2 + \mathbb{E}[X(s)] \mathbb{E}[X(t-s)] - \lambda^2 \int_0^s (1-F(s-r))dr \int_0^t (1-F(t-r))dr\\
\end{align*}
From here, we could use the knowledge that $X(s)$ and $X(t)$ are Poisson processes to write integrals for each term. I'm just not convinced that this approach is correct, namely because I don't see how to simplify it once I write down all the integrals. Am I on the right track? Any help would be fantastic. 
 A: Let $\{T_i:i=1,2,\ldots\}$ be the service times which are i.i.d. with common distribution function $F$ and $\{N(t):t\geqslant 0\}$ the arrival process. If customer $i$ enters the system at time $s\leqslant t$, then customer $i$ will remain in the system at time $t$ with probability
$$
p(s) = \mathbb P(T_i>t-s) = 1 - F(t-s).
$$
We require the following lemma:

Let $p:[0,t]\to[0,1]$ be an integrable function. Then $R(t)=\sum_{n=1}^\infty \mathsf 1_{(0,t]}(T_n)p(T_n)$ has Poisson distribution with parameter $\lambda\int_0^tp(s)\ \mathsf ds$.

Conditioned on $\{N(t)=m\}$ we have $R(t) = \sum_{i=1}^m I(T_i)$ where $I(T_i)$ has Bernoulli distribution with parameter $p(T_i)$. By Campbell's theorem, $R(t)\stackrel{\mathrm d}=\sum_{i=1}^mI(U_{(i)})$ where the $U_{(i)}$ are the order statistics of an i.i.d. sample $(U_1,\ldots,U_n)$ where $U_1\sim\mathrm{Unif}(0,t)$. Therefore $R(t)= \sum_{i=1}^m I(U_i)$, and the $I(U_i)$ are Bernoulli random variables with parameter $$\alpha = \frac 1t\int_0^t p(s)\ \mathsf ds. $$ Hence
$$
\mathbb P(R(t) = k\mid N(t)=m) = \binom mk \alpha^k (1-\alpha)^{m-k},
$$
so that
\begin{align}
\mathbb P(R(t) = k) &= \sum_{m=0}^\infty \mathbb P(R(t)=k\mid N(t)=m)\mathbb P(N(t)=m)\\
&= \sum_{m=k}^\infty \binom mk \alpha^k(1-\alpha)^{m-k} e^{-\lambda t}\frac{(\lambda t)^m}{m!}\\
&= e^{-\lambda t}\frac{\alpha^k}{k!} \sum_{m=k}^\infty (1-\alpha)^{m-k} \frac{(\lambda t)^m}{(m-k)!}\\
&= e^{-\lambda t}\frac{(\lambda\alpha t)^k}{k!} \sum_{m=0}^\infty \frac{(\lambda(1-\alpha)t)^m}{m!}\\
&= e^{-\lambda t}\frac{(\lambda\alpha t)^k}{k!} e^{\lambda(1-\alpha)t}\\
&= e^{-\lambda(1-\alpha)t}\frac{(\lambda\alpha t)^k}{k!},
\end{align}
so that $R(t)$ is Poisson distributed with parameter $\lambda\alpha t = \lambda\int_0^t p(s)\ \mathsf ds$.
Now, since $p(s) = 1-F(t-s)$ is a monotone bounded function, it is integrable, hence by the preceeding lemma $X(t)$ is a Poisson random variable with parameter
$$
\lambda\int_0^t p(s)\ \mathsf ds = \lambda \int_0^t (1-F(t-s))\ \mathsf ds = \lambda \int_0^t (1-F(s))\ \mathsf ds.
$$
Denote by $\lambda(t) = \lambda \int_0^t (1-F(s))\ \mathsf ds$. Then 
\begin{align}
\mathrm{Cov}(X(s),X(t)) &= \mathbb E[X(s)X(t)] - \mathbb E[X(s)]\mathbb E[X(t)]\\
&=\mathbb E[X(s)(X(t)-X(s)+X(s))] - \lambda(s)\lambda(t)\\
&= \mathbb E[X(s)^2] + \mathbb E[X(s)(X(t)-X(s))] - \lambda(s)\lambda(t)\\
&= \mathrm{Var}(X(s)) + \mathbb E[X(s)]^2 +\mathbb E[X(s)]\mathbb E[X(t)-X(s)] - \lambda(s)\lambda(t)\\
&= \lambda(s) + \lambda(s)^2 +\lambda(s)\lambda(t-s) -\lambda(s)\lambda(t)\\
&= \lambda(s)(1 + \lambda(s)+\lambda(t-s) -\lambda(t)).
\end{align}
