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}$.
Let's first consider $f(z)$ analytic inside and on a boundary $C$ of a simply connected region. Then Cauchy's integral formula is $f(z_0) = \frac{1}{2i\pi}\oint_C\frac{f(z)}{z-z_0}dz$.
Proof: Let $\Gamma$ be a circle of radius $\epsilon > 0$ inside $C$ with center $z_0$. Then $\oint_C\frac{f(z)}{z-z_0}dz = \oint_{\Gamma}\frac{f(z)}{z-z_0}dz$. The equation for a circle with center $z_0$ and radius $\epsilon$ is $\lvert z - z_0\rvert = \epsilon\iff z - z_0 = \epsilon e^{i\theta}$. Let $z = z_0 + \epsilon e^{i\theta}$ so $dz = i\epsilon e^{i\theta}d\theta$.
\begin{align}
\oint_{\Gamma}\frac{f(z)}{z-z_0}dz &= \int_0^{2\pi}\frac{f(z_0 + \epsilon e^{i\theta})}{\epsilon e^{i\theta}}i\epsilon e^{i\theta}d\theta\\
&= i\lim_{\epsilon\to 0}\int_0^{2\pi}f(z_0+\epsilon e^{i\theta})d\theta\\
&= 2i\pi f(z_0)
\end{align}
Thus, $f(z_0) = \frac{1}{2i\pi}\oint_C\frac{f(z)}{z-z_0}dz$.
Let $f(z)$ be analytic in a region bounded by two concentric circles $C_1$ and $C_2$ where the radius of $C_1$ is greater than the radius of $C_2$. By Cauchy's integral formula, we are able to show that $f(z_0) = \frac{1}{2i\pi}\oint_{C_1}\frac{f(z)}{z-z_0}dz-\frac{1}{2i\pi}\oint_{C_2}\frac{f(z)}{z-z_0}dz$.
Laurent's theorem: If $f(z)$ is analytic inside and on the boundary of an annular region bounded by two concentric circles centered at $z_0$ with radii $r_1$ and $r_2$, then for all $z$ in the annular region
$$
f(z) = \sum_{n=0}^{\infty}a_n(z-z_0)^n+\sum_{n=1}^{\infty}\frac{a_{-n}}{(z-z_0)^n}
$$
where the coefficients are defined as
\begin{align}
a_n &= \frac{1}{2i\pi}\oint_{C_1}\frac{f(w)}{(w-z)^{n+1}}dw\\
a_{-n} &= \frac{1}{2i\pi}\oint_{C_2}\frac{f(w)}{(w-z)^{-n+1}}dw
\end{align}
Proof: By Cauchy's integral formula, we have
$$
f(z) = \frac{1}{2i\pi}\oint_{C_1}\frac{f(w)}{w-z}dz-\frac{1}{2i\pi}\oint_{C_2}\frac{f(w)}{w-z}dz\tag{1}
$$
for $z$ in the annular region.
Since you are concerned with the $a_{-1}$, I am only going to show the work for the second integral in equation $(1)$. Consider $\frac{-1}{w-z}$.
\begin{align}
\frac{-1}{w-z} &= \frac{1}{z-z_0}\frac{1}{1 - \frac{w-z_0}{z-z_0}}\\
&= \frac{1}{z-z_0}\sum_{k=0}^{\infty}\Bigl(\frac{w-z_0}{z-z_0}\Bigr)^k\\
&= \frac{1}{z-z_0} + \frac{w-z_0}{(z-z_0)^2} + \frac{(w-z_0)^2}{(z-z_0)^3}+\cdots + \frac{(w-z_0)^{n-1}}{(z-z_0)^n}\\
&+ \text{higher order terms}\\
-\frac{1}{2i\pi}\oint_{C_2}\frac{f(w)}{w-z}dz &=
\frac{1}{2i\pi}\oint_{C_2}\frac{f(w)}{z-z_0}dw + \frac{1}{2i\pi}\oint_{C_2}\frac{w-z_0}{(z-z_0)^2}f(w)dw + \frac{1}{2i\pi}\oint_{C_2}\frac{(w-z_0)^2}{(z-z_0)^3}f(w)dw\\
&+\cdots + \frac{1}{2i\pi}\oint_{C_2}\frac{(w-z_0)^{n-1}}{(z-z_0)^n}f(w)dw + \text{higher order terms}\\
&= \frac{a_{-1}}{z-z_0} + \frac{a_{-2}}{(z-z_0)^2} + \frac{a_{-3}}{(z-z_0)^3}+\cdots + \frac{a_{-n}}{(z-z_0)^n}\\
&+ \text{higher order terms}
\end{align}
Then
$$
a_{-1}=\frac{1}{2i\pi}\oint_{C_2}f(w)dw.
$$
Best Answer
You are right about the value of the residue. And $0$ is an essential singularity indeed. Being an essential singularity means that, if your Laurent series is $\sum_{n=-\infty}^\infty a_nz^n$, then there are infinitely many negative integers $n$ with $a_n\ne0$. That's the case here; it's not zero for every negative odd integer.