Is the centered product of a Gaussian and Bernoulli r.v. sub-Gaussian?! Tailbound needed Let $X_1, X_2, \dots, X_n$ i.i.d. $N(0,\sigma^2)$ and let $Y_1, Y_2, \dots, Y_n$ be independent and identical Bernoulli random variables (where $Y_i$ may depend on $X_i$).
I am searching for a tailbound / concentration inequality of the form 
$$P(|\frac{1}{n}\sum_{i=1}^n X_iY_i  - \mathbb{E}(X_iY_i)|>z)\leq 2\exp(-c z^2)$$
for some specific value of $c$ (which should obviously depend on the variance of $X_iY_i$). 
If the tailbound holds I am very interested in a concrete (and "sharp") value of $c$!). 

Possible way to go: 
For me it looks like as if $X_iY_i-\mathbb{E}(X_iY_i)$ will still be sub-Gaussian (with which parameter!?). If this could be checked, then one would have to apply a Hoeffding bound to bound the sum of independent subgaussian random variables (which are sub-Gaussian again) and would be done.
However I am having trouble showing that $X_iY_i-\mathbb{E}(X_iY_i)$ is sub Gaussian and finding the correct parameter. 
Any help is greatly appreciated.

Edit: as Whuber pointed out with reference to wikipedia it is easy to see that $Z_i:=X_iY_i-\mathbb{E}(X_iY_i)$ is sub-Gaussian by checking the $\Psi_2$ condition. (Here done in more detail than needed, see Whuber's comment: it would have been enough to check the condition directly for $X_iY_i$)
Indeed: since $Z^2 \leq 2(X_iY_i)^2 + 2 \mathbb{E}(X_iY_i)^2$ and $|E(X_iY_i)|<d$ for some $d>0$ we have for all $a>0:$
\begin{align*}
\mathbb{E}(\exp(aZ^2)) & \leq \mathbb{E}(\exp(2a (X_iY_i)^2 + 2a\mathbb{E}(X_iY_i)^2))\\
& =\exp(2a\mathbb{E}(X_iY_i)^2)\mathbb{E}(\exp(2a (X_iY_i)^2))\\
& \leq \exp(2ad^2) \mathbb{E}(\exp(2a (X_iY_i)^2))\\
& \leq \exp(2ad^2) \mathbb{E}(\exp(2a X_i^2)) < \infty,
\end{align*}
since $X_i$ itself is sub-Gaussian and hence follows the $\Psi_2$ condition. 
Hence $Z_i:=X_iY_i -\mathbb{E}(X_iY_i)$ is sub-Gaussian with some parameter $b$ and $\sum_{i=1}^n Z_i$ is sub-Gaussian with parameter $nb$
However: I am still in need of a concrete value of the constant $c$ (or equivalently: $b$) (as sharp as possible) 
The true dependency between $X_i$ and $Y_i$ is too complicated to be given here, hence a more general bound on $c$ would be sufficient (which I feel should be possible, since |Y_i| is bounded by 1). 
However, if it is of any help/as a starter one could think that for each $i$, the relationship of $(X_i,Y_i)$ could be described by a logistic regression model. 
 A: For explicit absolute constants in the bound, you can use symmetrization and contraction. Indeed, if $\epsilon_1,...,\epsilon_n$ are independent random variables taking values $\pm1$ with probability $1/2$ (and independent of the $X_i$ and the $Y_i$, then by symmetrization (noting that $t\to\exp(\lambda t)$ is convex) we get
$E[ \exp(\lambda \sum_{i=1}^n (X_iY_i - E[X_iY_i]) \le E\exp(2\lambda \sum_{i=1}^n \epsilon_i X_i Y_i).$
Next, by contraction you get rid of the Bernoulli (even if $X_i,Y_i$ are dependent):
$E\exp(2\lambda \sum_{i=1}^n \epsilon_i X_i Y_i) \le E\exp(2\lambda \sum_{i=1}^n X_i)$.
But now $\sum_{i=1}^n X_i$ is $N(0, n\sigma^2)$ so the right hand side equals $\exp(2 \lambda^2 n \sigma^2)$. This implies by Markov inequality with $\lambda =t/(2\sigma\sqrt n)$
\[
P\left( \frac{1}{\sqrt n}\sum_{i=1}^n (X_i Y_i-E[X_i Y_i]) > t \sigma
\right) \le e^{-\lambda t \sigma \sqrt n + 2n\lambda^2\sigma^2 } = e^{-t^2/4}.
\]
We may use the same argument for the other direction.
What is symmetrization?
If $\epsilon_1,...,\epsilon_n$ are iid $\pm 1$ as above, independent of $Z_1,...,Z_n$ which are iid with $E[Z_i]=0$, then for any convex positive function $F$,
$E[ F(\sum_{i=1}^nZ_i) ]\le E[ F(\sum_{i=1}^n\epsilon_i Z_i) ]$.
Why is this true? of $Z_1',...,Z_n'$ are independent copies of $Z_1,...,Z_n$ then by Jensen's inequality
\begin{align}
E[ F(\sum_{i=1}^nZ_i) ] & \le E[ F(\sum_{i=1}^n (Z_i-Z_i') ] \\
&= E[ F(\sum_{i=1}^n \epsilon_i(Z_i-Z_i') ] \\
&\le \frac 1 2 E[ F(2\sum_{i=1}^n \epsilon_iZ_i) ] 
+ \frac 1 2 E[ F(2\sum_{i=1}^n -\epsilon_i Z_i') ] \\
&= E[ F(2\sum_{i=1}^n \epsilon_i Z_i) ].
\end{align}
Here, the first inequality is Jensen's and the $Z_i'$ have mean zero, the first equality is because $\epsilon_i(Z_i-Z_i')$ has the same distribution as $Z_i-Z_i'$ (because its distribution is symmetric), the second inequality is Jensen's again, and the last equality is because $\epsilon_iZ_i$ and $-\epsilon_iZ_i'$ are equal in distribution.
What is contraction?
If you have deterministic $x_1,...,x_n$ and random signs $\epsilon_1,...,\epsilon_n$ as above, then
\[
\sup_{\boldsymbol y\in [-1,1]^n} E_\epsilon [F(\sum_{i=1}^n \epsilon_i y_i x_i)
= E_\epsilon[F(\sum_{i=1}^n \epsilon_i x_i).
\]
(Here the only random variables are the $\epsilon_i$s).
Why is this true? Because the function $G:\boldsymbol y \to E[F(\sum_{i=1}^n \epsilon_i y_i x_i)$ from $[-1,1]^n\to R$ is convex, and convex functions over polytopes such as $[-1,1]^n$ attains their maximum at a vertex. Here, the vertices of $[-1,1]^n$ are of the form $y_i=\pm 1$ for each $i=1,...,n$. But for any $\boldsymbol y$ of such form, $G(\boldsymbol y) = E_\epsilon[F(\sum_{i=1}^n \epsilon_i x_i)$ because $(\epsilon_1,...,\epsilon_n)$ has the same distribution as $(y_1\epsilon_1,...,y_n\epsilon_n)$ when $y_i=\pm 1$.
Contraction for random $X_1,...,X_n$ and $Y_1,...,Y_n$...
Now if $X_1,...,X_n$ and $Y_1,...,Y_n$ are random with $Y_i\in[-1,1]$ (as in your example), then by conditioning on the $X_i,Y_i$'s,
\begin{align}
E[ F( 2 \sum_{i=1}^n \epsilon_i X_i Y_i) ]
& =
E\Big[ E[ F( 2 \sum_{i=1}^n \epsilon_i X_i Y_i) \Big|X_1,...,X_n,Y_1,...,Y_n] \Big] \\
&\le
E\Big[ \max_{\boldsymbol y\in[-1,1]^n} E[ F( 2 \sum_{i=1}^n \epsilon_i X_i y_i) \Big|X_1,...,X_n] \Big]  \\
& =
E\Big[ E[ F( 2 \sum_{i=1}^n \epsilon_i X_i) \Big|X_1,...,X_n] \Big]
= E[ F(2 \sum_{i=1}^n \epsilon_i X_i ) ].
\end{align}
A: I think I got one answer. It doubt it to be the sharpest bound, but at least it is a usable bound. 
I am using Theorem 2.1 from this rough draft: 
It is easy to see that (III) in Theorem 2.1 holds:
$$\mathbb{E}((X_iY_i-\mathbb{E}(X_iY_i))^{2k})\leq \frac{(2k)!}{2^k k!} (2\sigma)^{2k}=\frac{(2k)!}{2^k k!}  \theta^{2k},$$
where $\theta = 2\sigma$
From the proof $ (III) \rightarrow (I)$ given in Appendix A of the rough draft (p.39) we then may immediately conclude that for all $\lambda>0$ we have
$$\mathbb{E}(\exp(\lambda (X_iY_i-\mathbb{E}(X_iY_i))) \leq \exp(\frac{(\lambda \sqrt{2}\theta)^2}{2}),$$
i.e., by definition,  $Z_i:=(X_iY_i-\mathbb{E}(X_iY_i))$ is sub-Gaussian with parameter $\sqrt{2}\theta =  2^{\frac{3}{2}}\sigma$.
Since $Z_1,\dots Z_n$ are independent sub-Gaussian r.v. with parameter $2^{\frac{3}{2}}\sigma$. It follows from the Hoeffding bound that also $$P(\frac{1}{n}|\sum_{i=1}^n (X_iY_i-\mathbb{E}(X_iY_i))| \geq t) \leq 2 \exp(-\frac{t^2 n}{2^4\sigma^2})= 2\exp(-ct^2)$$
where $c \equiv c_n = \frac{n}{2^4 \sigma^2}$
(hence the value $c$ is not a pure constant, but depends also on $n$).
Moreover I think that one might still improve the factor $\frac{1}{2^4}$ appearing in $c$.
