Here, I'll prove that $\lim_{\epsilon\rightarrow 0} {\rm Im} f(\lambda - i \epsilon) = \lambda P(\lambda)$. No contour integration is involved. I'll also touch on why the imaginary part is only nonzero on the support of $P$, and why this is related to the pole in the integral defining $P$. The argument essentially involves proving this useful formula, which I'm assuming you don't know (that's why I prove it instead of just using it). At the very end, I'll also prove the other part of the claim, relating $f(z)$ to $\sum_n R_n z^{-n}$.
Let's first expand the right hand side using $P=\sum_i \delta(\lambda-\lambda_i)$
\begin{equation}
f(z) = \frac{1}{\pi} \int d \lambda \frac{\lambda P(\lambda)}{z-\lambda} = \frac{1}{\pi} \sum_i \int d \lambda \frac{\lambda \delta(\lambda_i - \lambda)}{z-\lambda}= \frac{1}{\pi}\sum_i \frac{\lambda_i}{z - \lambda_i}
\end{equation}
Now, let's write $z=x - i \epsilon$. For now, $x$ is just some arbitrary real number, and $\epsilon$ is a real number that we will send to be zero. Then we can rewrite $f(z)$ as
\begin{eqnarray}
f(x-i\epsilon) &=& \frac{1}{\pi} \sum_i \frac{\lambda_i}{(x - \lambda_i) - i \epsilon} \\
&=& \frac{1}{\pi} \sum_i \frac{\lambda_i(x-\lambda_i)}{(x - \lambda_i)^2 + \epsilon^2} + \frac{i}{\pi} \sum_i \frac{\epsilon \lambda_i }{(x-\lambda_i)^2 + \epsilon^2}
\end{eqnarray}
and therefore
\begin{equation}
{\rm Im} f(z) = \frac{1}{\pi} \sum_i \frac{\epsilon \lambda_i}{(x-\lambda_i)^2 + \epsilon^2}
\end{equation}
Now we take the limit $\epsilon\rightarrow 0$. If $x\neq \lambda_i$, then clearly $\lim_{\epsilon \rightarrow 0} {\rm Im} f(x - i \epsilon) = 0$.
However, if $x = \lambda_j$, $\lim_{\epsilon\rightarrow 0} {\rm Im} f(x-i\epsilon)=0/0$, an indeterminate form, so we need to be more careful. First, let's get back control by moving slightly away from the dangerous point $\lambda_j$ by writing $x = \lambda_j + y$ (we'll set $y=0$ at the end). Then, note that every term in the sum except $i=j$ vanishes, by the argument above (here I'm assuming the spectrum is non-degenerate; I'll leave it as an exercise to generalize this to the case with degenerate eigenvalues :-)). As a result, we have
\begin{eqnarray}
\lim_{\epsilon \rightarrow 0} {\rm Im} f(\lambda_j - i \epsilon) &=& \lim_{y\rightarrow 0} \lim_{\epsilon \rightarrow 0} {\rm Im} f(\lambda_j + y - i \epsilon) \\
&=& \lim_{y \rightarrow 0} \lambda_j \times \left( \lim_{\epsilon \rightarrow 0} \frac{1}{\pi} \frac{\epsilon}{y^2 + \epsilon^2} \right) \\
&=& \lim_{y \rightarrow 0} \lambda_j \delta(y) \\
&=& \lambda_j \delta(0)
\end{eqnarray}
where we have used one of the representations of the delta function (see "Poisson kernel" on wikipedia). In typical physicist fashion I've been a bit cavalier about the order of limits, although it works out here.
Unraveling what we've done, we can now say that
\begin{equation}
\lim_{\epsilon \rightarrow 0} {\rm Im} f(z - i \epsilon) =
\cases{
\lambda_j \delta(0) = \lambda_j P(\lambda_j), \ \ & $z=\lambda_j$ \\
0, \ \ & otherwise \\
}
\end{equation}
where I used $P(\lambda_j) = \delta(0)$, which follows immediately from the definition of $P(\lambda)$ and the fact that $\lambda_j$ is an eigenvalue.
In other words, we have shown that the imaginary part of $f$ is non-zero only on the support of $P$, as desired. This proves
\begin{equation}
\lim_{\epsilon \rightarrow 0} {\rm Im} f(\lambda-i\epsilon) = \lambda P(\lambda).
\end{equation}
By the way, as a byproduct of our argument, we've essentially proven the following curious but useful formula
\begin{equation}
\lim_{\epsilon \rightarrow 0} \frac{1}{z - i \epsilon} = \mathcal{P} \frac{1}{z} + i \pi \delta(z)
\end{equation}
where $\mathcal{P} \frac{1}{z}$ is the principal part of $1/z$. If we had started from this formula, we could have done the argument in one line, as follows
\begin{equation}
\lim_{\epsilon \rightarrow 0} {\rm Im} f(z - i \epsilon) = \int d \lambda \lambda P(\lambda) \left[ {\rm Im} \lim_{\epsilon \rightarrow 0} \frac{1}{\pi} \frac{1}{z-\lambda - i\epsilon} \right] = \int d \lambda \lambda P(\lambda) \delta(z-\lambda) = zP(z).
\end{equation}
This formula also makes it clear that the imaginary part of $f$ comes from the pole, which was one of the key points of your question.
There's also another part of the claim, namely that
\begin{equation}
\sum_{n=1}^\infty R_n z^{-n} = \int d \lambda \frac{\lambda P(\lambda)}{z-\lambda}
\end{equation}
where
\begin{equation}
R_n = \sum_i \lambda_i^n
\end{equation}
We can prove this straightforwardly by starting from the left hand side
\begin{eqnarray}
\sum_{n=1}^\infty R_n z^{-n} &=& \sum_{n=1}^\infty \sum_i \lambda_i^n z^{-n} = \sum_i \left[ \sum_{n=0}^\infty \left(\frac{\lambda_i}{z}\right)^n - 1 \right]= \sum_i \left[\frac{1}{1 - \lambda_i/z} - 1\right]\\
& =& \sum_i \frac{\lambda_i}{z - \lambda_i} \\
&=& \int d \lambda \frac{\lambda P(\lambda)}{z-\lambda}
\end{eqnarray}
where we used $P(\lambda) = \sum_i \delta(\lambda-\lambda_i)$. We've had to assume that $|\lambda_i/z| < 1$ here.
