Let $\psi(x) = \exp(-x^2)$. The (general) Wavelet transform $W: L^2(\mathbb R) \rightarrow X$ is defined by $$Wf(a, b) = b^{-1/2}\int_{x\in \mathbb R} f(x) \psi\left( \frac{x-a}{b}\right)dx$$ for $(a, b) \in \mathbb R_+ \times \mathbb R$.
My question is: What is the functional space $X$?
NOTE:
Note that the function $\psi$ above is a Gaussian-like function, which is not exactly a wavelet function. In general, a function $\psi$ is called a wavelet function if it satisfies 1) $\psi$ has compact support or at least fast decaying at infinity and 2) $\int_{x \in \mathbb R} \psi(x) dx=0$. The given function $\psi$ however only satisfies the first requirement. We therefore refer it as a "general" Wavelet transformation.
To have a clue for our question. Let's us consider a particular case, when $b=b_0$ is constant, then $Wf$ is just a function of $a$. In this case $$Wf \in X = C_0(\mathbb R).$$ Indeed, in this case, one can see that $Wf = b_0^{-1/2} f\star \psi(\cdot/b_0)$, i.e. $W$ is a simply a convolution operator. A (well-known?) fundamental result states that: The convolution operator of two $L^2$ functions is $C_0$. Hence $ Wf\in C_0(\mathbb R)$.
For the original question, numerical simulations also suggest that $X=C_0(\mathbb R_+\times \mathbb R)$. It it correct?
Since $X = C_0(\mathbb R)$ satisfies the letter of your question, I assume you want the smallest such $X$, or in other words, the range of $W$.
First, your reasoning also applies when considering $Wf$ as a bivariate function, so you can in fact prove $Wf \in C_0(\mathbb R \times (0,\infty))$ (you have to exclude the case $b=0$).
There is more though: because your kernel is $C^\infty$ and all its derivatives are $L^\infty$, you can in fact show that $Wf$ is $C^\infty$. More than that, it is holomorphic on $\mathbb C \times \{z \in \mathbb C : \mathrm{Re}(z) > 0 \}$, see below.
Of course not all holomorphic functions on this domain are in the range of $Wf$; for example, by the Cauchy-Schwarz inequality $$ |Wf(a,b)| \leq \left(\int |f|^2 \right)^{1/2} $$ for every $(a,b)\in \mathbb R \times (0,\infty)$, so your function $Wf$ must be bounded on $\mathbb R \times (0,\infty)$ (it cannot be bounded as a complex function though). To establish the range of $W$, define an inverse transform (that recovers $f$ from $Wf$) and see which analytic functions are sent to $L^2$ functions by this inverse transform.
As a final note, there is always the possibility that $L^2$ is not the best space to consider; for example, the Fourier transform is defined on tempered distributions, which are a much larger class of objects.
Proof of holomorphicity (you can skip)
For simplicity I will only prove it when $a$ and $b$ are both real. To prove that a function is holomorphic, we need to prove that it is equal to its infinite Taylor series near any point $(a,b)$ of its domain of definition.
The first thing I want to do is define $Vf(a,b) = Wf(a,b^{-1})$. This makes the computations a bit easier without changing the conclusion.
Fix $a\in\mathbb R$ and $b>0$, and reparametrize it by $H(u,v) = Vf(a+u, b+v)$. Then \begin{align} H(u,v) &= (b+v)^{-1/2} \int f(x) e^{-(b+v)^2 (x-a-u)^2} dx \\ &= (b+v)^{-1/2} e^{-(b+v)^2 (a+u)^2}\int f(x) e^{-(b+v)^2 \left[ x^2 - (a+u)x \right]} dx \end{align} The function $v \mapsto (b+v)^{-1/2} e^{-(b+v)^2 (a+u)^2}$ is indeed holomorphic for $|v|<b$. If we prove that the integral is holomorphic as a function of $(u,v)$, the conclusion follows. This integral equals \begin{align} &\int f(x) e^{-b^2 x^2} e^{-((b+v)^2+b^2)x^2 + (b+v)^2 (a+u)x} dx \\ &= \int \sum_{n\geq 0} f(x) e^{-b^2 x^2} \frac 1 {n!} \left(\left(-(b+v)^2+ b^2\right)x^2 + (b+v)^2 (a+u)x\right)^n dx \\ &= \int \sum_{n,j,k\geq 0} f(x) e^{-b^2 x^2} c_{n,j,k} u^j v^k P_{n,j,k}(x) dx . \qquad (1) \end{align} where the $P_{n,j,k}$ are polynomials in $x$: the first line is the expansion of the exponential, the second is the binomial theorem – and by writing down the binomial theorem explicitely, we can prove the very useful fact that for every $u,v,x$ with $|v|<b$, \begin{align} & \sum_{j,k\geq 0} \left| c_{n,j,k} u^j v^k P_{n,j,k}(x) \right| \leq \frac 1 {n!} \left(\left(b^2-(b-|v|)^2\right)x^2 + (b+|v|)^2 (|a|+|u|)|x|\right)^n . \end{align} Can we exchange the sums and integral in $(1)$? Let us check the condition of Fubini’s theorem: for $|v|<b$, \begin{align} &\int \sum_{j,k,n\geq 0} |f(x)| e^{-b^2 x^2} \left| c_{n,j,k} u^j v^k P_{n,j,k}(x) \right| dx \\ &\leq \int \sum_{n\geq 0} |f(x)| e^{-b^2 x^2} \frac 1 {n!} \left(\left(-(b+|v|)^2+ b^2\right)x^2 + (b+|v|)^2 (|a|+|u|)|x|\right) dx \\ &= \int |f(x)| \exp\left( -b^2 x^2+ |v|(2b+|v|) x^2 + (b+|v|)^2 (|a|+|u|)|x| \right) dx \\ &\leq \left( \int |f|^2 \right)^{1/2} \left(\int \exp\left( -2(b-|v|)^2 x^2+ 2 (b+|v|)^2 (|a|+|u|)|x| \right) dx \right)^{1/2} \end{align} where we used the Cauchy-Schwarz inequality in the last step. We can check that the second integral is finite, meaning the function in $(*)$ is integrable, so we can exchange the sums and integral in any way we want: $$ (1) = \sum_{j,k\geq 0} h_{j,k} u^j v^k $$ with $$ h_{j,k} := \sum_{n\geq 0}\int f(x) e^{-b^2 x^2} c_{n,j,k} P_{n,j,k}(x) dx . $$ This means that $(1)$, and thus $H$, is holomorphic for $(u,v)$ close to zero, or that $Vf$ is holomorphic close to $(a,b)$. Since this is true for every $(a,b) \in \mathbb R \times (0,\infty)$, it means that $Vf$ is holomorphic on a neighborhood of that set.