The Laurent series in any given annulus is uniquely determined. (You should look up a proof of this. It's important.) In this case, you have convergence in a punctured disk. To determine this disk, look at the singularities of the function and note the distance from 1 to the nearest singularity. (If you return to the proof of the existence of Laurent series, you will see that the radius of convergence will be the largest possible punctured disk around with point that doesn't contain another singularity.)
The hard part about finding the Laurent series here is inverting the denominator. If you write out the power series for the denominator, you will see that it is
$$(z-1)^2(a+b(z-1)+c(z-1)^2+\dots)$$
where $a,b,c$ are some constants. So you can factor out the $\frac{1}{(z-1)^2}$ term and use the normal procedure for inverting power series with a nonzero constant term. You can use your idea (writing everything in terms of $z-1$) to take care of the numerator.
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}$.
Best Answer
Assume wlog $z_0=0$ and note that $f(z)=\sum_{n >0}b_nz^{-n}+g(z)$ where $g$ is holomorphic on a disc of radius $\delta$ around $0$, and let $h(z)=\sum_{n >0}b_nz^{-n}$.
For any $r< \delta$ let $\max_{|z|=r} |f(z)|=M_r$. Then for any $n \ge 1$ we have $$|b_n|=\frac{1}{2\pi i}\int_{|z|=r}f(z)z^{n-1}dz$$ and since $|dz|=r$ we have that $|b_n| \le M_r r^n$ hence $|b_n|^{1/n} \le M_r^{1/n}r$ so $\limsup |b_n|^{1/n} \le r$
Since $r>0$ can be taken as small as we want, we get that $\limsup |b_n|^{1/n} =0$ so $h(1/z)=\sum_{n \ge 1}b_nz^n$ is entire, which means precisely that $h$ is analytic on $\mathbb C\setminus \{z_0\}$ and we are done!