Here is another proof using complex analysis that avoids using a Laurent series expansion. However, this argument is not as direct as the argument using Laurent series expansions in the answer by Just a user.
We assume $A \neq 0$. Let $f\in \mathcal{L}(X)^{*}$. Let $\varepsilon > 0$. Note that if $\lambda\in\mathbb{C}\setminus\{0\}$ and $|\lambda | < (r(A) + \varepsilon )^{-1}$ then $r(A) + \varepsilon < |\lambda^{-1}|$, which implies that $I - \lambda A = \lambda (\lambda^{-1} I - A)$ has a bounded inverse. Hence we define $\Omega := \{\lambda\in\mathbb{C} : |\lambda | < (r(A) + \varepsilon )^{-1} \}$ and $u\colon \Omega\to\mathbb{C}$ by
\begin{equation}
u(\lambda) := \lambda f( (I - \lambda A)^{-1}) .
\end{equation}
The map $\lambda \mapsto f((\lambda^{-1}I - A)^{-1})$ is holomorphic on $\Omega\setminus\{0\}$ because it is the composition of the holomorphic map $\lambda\mapsto \lambda^{-1}$ on $\Omega$ with the holomorphic map $\lambda\mapsto f((\lambda I - A)^{-1})$ on $\{\lambda\in\mathbb{C} : |\lambda | > r(A) + \varepsilon \}$, so we see that $u$ is holomorphic on $\Omega\setminus\{0\}$. But also for each $\lambda\in\mathbb{C}$ with $|\lambda | < \|A\|^{-1}$ we have $\|\lambda A\| < 1$ and hence
\begin{equation}
u(\lambda ) = \lambda f \Big(\sum_{k=0}^{\infty} \lambda^{k} A^{k} \Big) = \sum_{k=0}^{\infty} \lambda^{k+1} f(A^{k}) . \tag{1}
\end{equation}
It follows that $u$ is holomorphic on $\{\lambda\in\mathbb{C} : |\lambda | < \|A\|^{-1}\} \cup (\Omega \setminus \{0\}) = \Omega$. Hence $u$ has a unique power series expansion on the open disk $\Omega$. Since the coefficients of the power series expansion of $u$ on $\Omega$ are determined by a power series expansion of $u$ on any smaller open disk centred at $0$, we see from $(1)$ that we have
\begin{equation}
u(\lambda ) = \sum_{k=0}^{\infty} \lambda^{k+1} f(A^{k}) \tag{2}
\end{equation}
for all $\lambda\in\Omega$.
Now let $\lambda\in\mathbb{C}$ with $|\lambda | > r(A) + \varepsilon$. Then $|\lambda^{-1}| < (r(A) + \varepsilon )^{-1}$, so $\lambda^{-1} \in \Omega$ and by $(2)$ we have
\begin{equation}
f((\lambda I - A)^{-1} ) = \lambda^{-1} f((I - \lambda^{-1} A)^{-1} ) = u(\lambda^{-1} ) = \sum_{k=0}^{\infty} \lambda^{-(k+1)} f(A^{k}) .
\end{equation}
As $\varepsilon > 0$ was arbitrary we obtain
\begin{equation}
f((\lambda I - A)^{-1} ) = \sum_{k=0}^{\infty} \lambda^{-(k+1)} f(A^{k}) \tag{3}
\end{equation}
for all $\lambda\in\mathbb{C}$ with $|\lambda | > r(A)$. This shows the first identity.
We now show the second identity. Let $f\in\mathcal{L}(X)^{*}$ and $\varepsilon > 0$. Let $\alpha\in\mathbb{R}$ with $r(A) < \alpha < r(A) + \varepsilon$. By $(3)$ the series $(\sum_{k=0}^{n} \alpha^{-(k+1)} f(A^{k}))_{n\in\mathbb{N}}$ converges. This implies that the sequence $(\alpha^{-(k+1)} f(A^{k}))_{k\in\mathbb{N}}$ is bounded. Hence there is some $M\in (0, \infty )$ such that
\begin{equation}
|(r(A) + \varepsilon )^{-(k+1)} f(A^{k})| \leq M \Big( \frac{\alpha}{r(A) + \varepsilon} \Big)^{k+1}
\end{equation}
for all $k\in\mathbb{N}$ and consequently the series
\begin{equation}
\Big( \sum_{k=0}^{n} (r(A) + \varepsilon )^{-(k+1)} f(A^{k}) \Big)_{n\in\mathbb{N}}
\end{equation}
converges absolutely. This shows the second identity.
It is also possible to show the first identity by using Gelfand's formula directly combined with this result as done in the answer by Yiorgos S. Smyrlis. The potential issue is that the argument is circular if you are trying to use the result in the question to show Gelfand's formula. In particular, you cannot use Gelfand's formula if you are using the result in the question to show $\limsup_{k\to\infty} \|A^{k}\|^{\tfrac{1}{k}} \leq r(A)$. However, there are other ways to derive Gelfand's formula that avoid the result in the question. Such proofs are mentioned in the answers to this question. One proof involves writing powers of $A$ in terms of a contour integral using the Cauchy integral formula or holomorphic functional calculus and then estimating this integral. Proofs of this argument can be found in Theorem 11.4 of Linear Functional Analysis by Alt or Chapter 6 of Functional Analysis by van Neerven. There is also a proof that avoids complex analysis entirely which can be found in the discussion following Theorem 1.2.8 of A Course in Commutative Banach Algebras by Kaniuth. To sketch the argument, a short argument demonstrates it is sufficient to show that the assumptions $\lim_{k\to\infty} \|A^{k}\|^{\tfrac{1}{k}} = 1$ and that $\lambda I - A$ has a bounded inverse for every $\lambda\in\mathbb{C}$ with $|\lambda | \geq 1$ lead to a contradiction. It is now shown for every $n\in\mathbb{N}$ that $I - A^{n}$ has a bounded inverse. This step makes use of existence of $n$-th roots of unity which do not exist over the real field in general. There is then an argument around a page long using elementary results to show that $(A^{k})_{k\in\mathbb{N}}$ converges to $0$ which contradicts that $\lim_{k\to\infty} \|A^{k}\|^{\tfrac{1}{k}} = 1$ and the result follows. These proofs just mentioned can also be generalised to the setting of Banach algebras.
Note that the last proof discussed above makes use of complex numbers in an essential way, since Gelfand's formula and the result in question fail if $X$ is a Banach space over the real field. To see why, note that it is possible for the spectrum of a bounded linear operator on $X$ to have empty spectrum. However, even when the spectrum is nonempty the results can still fail. Take $X = \mathbb{R}^{3}$ and $B\colon \mathbb{R}^{3}\to\mathbb{R}^{3}$ defined by $B(x_{1}, x_{2}, x_{3}) := (0, -x_{3}, x_{2})$. Then $\sigma (B) = \{0\}$ while $\|B^{k}\| = 1$ for each $k\in\mathbb{N}$, so it follows that
\begin{equation}
r(B) = 0 \neq 1 = \lim_{k\to\infty} \|T^{k}\|^{\tfrac{1}{k}} .
\end{equation}
Furthermore, we see that $I-B$ has an inverse while the series $(\sum_{k=0}^{n} B^{k} )_{n\in\mathbb{N}}$ does not converge.
