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) \\ &=(2\pi)^{-k/2}\int_0^\infty\exp\big(c\,(\tfrac1{2\sqrt k})^kr^k-r^2/2\big)\,dr\\ &\times P\big(U_i>\tfrac1{2\sqrt k}\ \forall i\in[k]\big) =\infty\end{aligned}$$$$\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$.