You can obtain the desired result using absolute convergence and some observations about the annulus of convergence of Laurent series.
If $\{x_n\}_{n \in \mathbb{Z}}$, then let $R_+(\{x_n\}) = ( \limsup_{n \to \infty} \sqrt[n]{|x_n|} )^{-1}$, and $R_-(\{x_n\}) = \limsup_{n \to \infty} \sqrt[n]{|x_{-n}|}$.
Two relevant results (all summations are over $\mathbb{Z}$):
(i) If $\sum_m \sum_n |a_{m,n}| < \infty$, then the summation can be rearranged. In particular, $\sum_m \sum_n a_{m,n} = \sum_k \sum_l a_{l,k-l}$. (Since $\phi(m,n) = (m,m+n)$ is a bijection of $\mathbb{Z}^2$ to $\mathbb{Z}^2$.)
(ii) If $\sum_{n} |x_n| < \infty$, then $R_+(\{x_n\}) \ge 1$. Similarly, $R_-(\{x_n\}) \le 1$.
Let $f(z) = \sum_n f_n (z-z_0)^n$, $g(z) = \sum_n g_n (z-z_0)^n$. Let $R_+ = \min(R_+(\{f_n\}),R_+(\{g_n\}))$, $R_- = \max(R_-(\{f_n\}),R_-(\{g_n\}))$.
Choose $R_-<r < R_+$. Then $\sum_m |f_m| r^m$ and
$\sum_n |g_n| r^n$ are absolutely convergent sequences, and so $\sum_m \sum_n |f_m||g_n| r^{m+n} < \infty$. From (i) we have $\sum_m \sum_n |f_m||g_n| r^{m+n} = \sum_k (\sum_l |f_l| |g_{k-l}|) r^k < \infty$.
If we let $c_k = \sum_l f_l g_{k-l}$, this gives $\sum_k |c_k| r^k < \infty$, and so
$R_+(\{c_kr^k\}) = \frac{1}{r} R_+(\{c_k \}) \ge 1$, or, in other words, $R_+(\{c_k \}) \ge r$. Since $r<R_+$ was arbitrary, it follows that $R_+(\{c_k \}) \ge R_+$. The same line of argument gives $R_-(\{c_k \}) \le R_-$.
Hence we have $c(z) = f(z)g(z) = \sum_n c_n (z-z_0)^k$ on $R_- < |z-z_0| < R_+$, where $c_k = \sum_l f_l g_{k-l}$.
I wanted to prove the convergence for both parts of the series, because if both parts converge, then also the whole Laurent series.
This is a very good idea, indeed:
The power series $\sum_{n=0}^\infty a_n(z-z_0)^n$ converges in the disk $|z-z_0|<R$ where $R$ is given by the Cauchy-Hadamard formula $$R^{-1}= \limsup |a_n|^{1/n} $$
The series of negative terms can be transformed into a power series with positive powers using the substitution $\frac{1}{z-z_0}=w$. The series then becomes $$\sum_{n=1}^\infty a_{-n} w^n $$ and it converges in the disk $|w|<r$, where, again by Cauchy-Hadamard's formula $r^{-1}=\limsup |a_{-n}|^{1/n}$. Clearly the disk $|w|<r$ is equivalent to the annulus $\left| \frac{1}{z-z_0}\right|<r$.
Combining the above results, the problem is solved.
Best Answer
I'm going to take $z_0=0$ for simplicity. For $\rho_1, \rho_2$ with $r_1 < \rho_1 < \rho_2 < r_2$, and for $z$ with $\rho_1 < |z| < \rho_2$, Cauchy's integral formula applied to the annulus $\{ \rho_1 \leq |w| \leq \rho_2\}$ gives us that
$$f(z)= \frac{1}{2\pi i} \int_{|w|=\rho_2} \frac{f(w)}{w-z} \ \mathrm{d}w - \frac{1}{2 \pi i} \int_{|w|=\rho_1} \frac{f(w)}{w-z} \ \mathrm{d}w.$$
In the first integral we can write
$$\frac{f(w)}{w-z} = \sum_{n=0}^\infty \frac{z^nf(w)}{w^{n+1}},$$
and for each $z$ with $|z| < \rho_2$ this sum converges aboslutely uniformly on the circle $\{|w|=\rho_2\}$.
[Added in response to comment: there exists a positive real number $M$ such that $|f(w)| \leq M$ on $\{|w|=\rho_2\}$. Then for all $w, z$ with $|w|=\rho_2$ and $|z|<\rho_2$, and all $n \geq 0$, we have
$$\Bigg|\frac{z^nf(w)}{w^{n+1}}\Bigg| \leq \frac{M}{\rho_2} \Bigg(\frac{|z|}{\rho_2}\Bigg)^n.$$
The right-hand side converges (it's a geometric series with ratio $|z|/\rho_2 <1$) so the claimed sum converges absolutely uniformly on the circle by the Weierstrass M-test.]
So the first integral is equal to
$$\sum_{n=0}^\infty z^n \frac{1}{2\pi i} \int_{|w|=\rho_2} \frac{f(w)}{w^{n+1}} \ \mathrm{d}w,$$
and this converges and is valid on the whole of $\{|z| < \rho_2\}$. Letting $\rho_2 \rightarrow r_2$, we see that the positive part of the Laurent series converges on $\{|z| < r_2\}$.
For the negative part, do a similar expansion with $w^n/z^{n+1}$.