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}$.
The function $f(z)=\exp(z+1/z)$ is analytic in the annular region
$0<|z|<\infty$.
We can represent the factors $\exp(z)$ and $\exp(1/z)$ as Laurent
series over the annulus $0<|z|<\infty$, viz.,
$\exp(z)=\sum_{k=-\infty}^\infty a_k z^k$ and
$\exp(1/z)=\sum_{k=-\infty}^\infty b_k z^k$, where
$$a_k=\begin{cases} 1/k!&\text{for $k\ge0$}\\
0&\text{for $k<0$}\,,\end{cases}$$ and
$$b_k=\begin{cases} 0&\text{for $k>0$}\\
1/(-k)!&\text{for $k\le0$}\,.\end{cases}$$
The Laurent series of the product is
$f(z)=\sum_{k=-\infty}^\infty c_k z^k$, where
$c_n=\sum_{k=-\infty}^\infty a_k b_{n-k}$.
The latter simplifies to $c_n=\sum_{k=0}^\infty\frac1{k!} \frac1{(n+k)!}$.
Further, by comparison with the series expression for the modified Bessel
function of the first kind,
$$I_n(z)=\sum_{k=0}^\infty{(\frac12 z)^{n+2k}\over k!\,\Gamma(n+k+1)}\,,$$
it follows that $c_n=I_n(2)$.
Thus we have that
$$\exp\left(z+\frac1z\right)=
I_0(2)+\sum_{n=1}^\infty I_n(2)\left(z^n+{1\over z^n}\right).$$
Best Answer
First, you can commute the two sums to disentangle their indices in the following way : $$ \sum_{k\ge0}\sum_{0\le l\le k} \frac{(-1)^l}{(k-l)!} z^{n-k+2l} = \sum_{l\ge0}\sum_{k\ge l} \frac{(-1)^l}{(k-l)!} z^{n-k+2l} = \sum_{l\ge0}\sum_{k\ge0} \frac{(-1)^l}{k!} z^{n-k+l} $$ Now, the coefficient $a_{-1}$ will be associated to the case $n-k+l = -1$, hence $k = n+l+1$ and $$ \begin{array}{rcl} a_{-1} &=& \displaystyle \sum_{l\ge0} \frac{(-1)^l}{(n+l+1)!} \\ &=& \displaystyle \sum_{0\le l\le n} \frac{(-1)^l}{(n+l+1)!} + \sum_{l\ge n+1} \frac{(-1)^l}{(n+l+1)!} \\ &=& \displaystyle \sum_{0\le l\le n} \frac{(-1)^l}{(n+l+1)!} + \sum_{l\ge0} \frac{(-1)^{n+l+1}}{l!} \\ &=& \displaystyle \frac{(-1)^{n+1}}{e} + \sum_{0\le l\le n} \frac{(-1)^l}{(n+l+1)!} \end{array} $$ Unfortunately, the partial sum doesn't seem to possess a closed form.