This sub-Gaussian technique cannot work in general for any $k\ge3$. In particular, if $k\ge3$ and $X_1,\dots,X_k$ are iid standard normal, then 
$$E\exp\Big(c\prod_1^k X_i)=\infty$$
for any $c>0$. 

Indeed, let $[k]:=\{1,\dots,k\}$, $X:=(X_1,\dots,X_k)$, let $(U_1,\dots,U_k)$ denote a uniformly distributed unit random vector, and let $|\cdot|$ denote the Euclidean norm. Then for some real $C_k>0$
$$\begin{aligned}E\exp\Big(c\prod_1^k X_i\Big)
&\ge E\exp\Big(c\prod_1^k X_i\Big)1\big(X_i>\tfrac1{2\sqrt k}\,|X|\ \forall i\in[k]\big) \\ 
&=C_k\int_0^\infty\exp\big(c\,(\tfrac1{2\sqrt k})^kr^k-r^2/2\big)\,r^{k-1}dr\\
&\times P\big(U_i>\tfrac1{2\sqrt k}\ \forall i\in[k]\big)
=\infty\end{aligned}$$
if $c>0$ and $k\ge3$.