Moment matching: construction of a mixture of Gaussian distribution with lower moments identical to Gaussian This is a question related to the statistical model behind independent component analysis (ICA). 
We assume that $Z \sim N(0,1)$. Our goal is to construct a random variable $X$ that follows a Gaussian mixture distribution and satisfies 
$$
E[X^j] = E[Z^j]
$$ 
for $j \leq k-1$, while 
$$
E[X^k] - E[Z^k] = \epsilon.
$$ 
Here $k \geq 2$ is an integer.
My question is: Is there a principled way to construct the Gaussian mixture random variable $X$?
For $k = 4$, one possible construction is $X\sim 1/2 \cdot N(0,1+\varepsilon^{1/2}/\sqrt{3}) + 1/2 \cdot N(0,1-\varepsilon^{1/2}/\sqrt{3})$. However, is there a construction for general $k$?
 A: If non-zero means are allowed, the answer is yes (and it is based on quite different considerations, and so, is presented separately from the previous answer). 
Indeed, let $Z$ be a standard normal random variable (r.v.). 
Let $P_k$ denote the set of all probability measures $\nu$ on $\mathbb R$ such that $\int_{\mathbb R}x^j\nu(dx)=EZ^j$ for all $j=0,\dots,k$. Let $M_k:=\sup\{\int_{\mathbb R}x^{k+1}\nu(dx)\colon\nu\in P_k\}$ and 
$m_k:=\inf\{\int_{\mathbb R}x^{k+1}\nu(dx)\colon\nu\in P_k\}$. Later it will be shown that 
$$m_k<EZ^{k+1}<M_k,\tag{-1}
$$ 
with an emphasis on the strictness of the inequalities. 
Take any 
$$c\in(m_k,M_k). \tag{0} $$ 
Let us then show that there is a mixture $\mu_c$ of normal distributions such that $\int_{\mathbb R}x^{k+1}\mu_c(dx)=c$ and $\mu_c\in P_k$, that is, $\int_{\mathbb R}x^j\mu_c(dx)=EZ^j$ for all $j=0,\dots,k$. 
For $t\in(0,1)$, let 
$$c_t:=c-(1- t^{(k+1)/2})EZ^{k+1}. 
$$
Since $c_t\to c$ as $t\uparrow1$, one can choose $t$ close enough to $1$ so that, by $(0)$, 
$$c_t\in(t^{(k+1)/2}m_k,t^{(k+1)/2}M_k).$$
So, by the convexity of the set $P_k$ and my answer to the question at 
[mathoverflow.net/questions/229646]
there is a r.v. $X_t$ taking only finitely many values such that 
$$EX_t^j=E(\sqrt tZ)^j\text{ for }j=0,\dots,k\text{ and }EX_t^{k+1}=c_t. \tag{1}
$$ Without loss of generality, $X_t$ is independent of $Z$. Let $Y_t:=X_t+\sqrt{1-t}Z$. Then the distribution of $Y_t$ is a mixture of normal distributions (each of those distributions with variance $1-t$). 
Moreover, for $j=0,1,\dots$ 
$$EY_t^j=\sum_{i=0}^j\binom ji EX_t^i\, E(\sqrt{1-t}Z)^{j-i}. \tag{2} 
$$
On the other hand, $Z$ is equal in distribution to $\sqrt t Z_1+\sqrt{1-t}Z_2$, where $Z_1$ and $Z_2$ are independent standard normal r.v.'s. So, 
$EZ^j=\sum_{i=0}^j\binom ji E(\sqrt t Z)^i\, E(\sqrt{1-t}Z)^{j-i}$ 
for all $j=0,1,\dots$. 
Comparing this with $(2)$ and recalling $(1)$, we see that $EY_t^j=EZ^j$ for $j=0,\dots,k$, whereas 
$$EY_t^{k+1}=EX_t^{k+1}+\sum_{i=0}^k\binom{k+1}i EX_t^i\, E(\sqrt{1-t}Z)^{k+1-i}$$
$$=EX_t^{k+1}+\sum_{i=0}^k\binom{k+1}i E(\sqrt t Z)^i\, E(\sqrt{1-t}Z)^{k+1-i}$$
$$
=c_t+EZ^{k+1}-E(\sqrt t Z)^{k+1}=c, 
$$
as desired.
Addendum. 
As promised, I am adding the proof of inequalities $(-1)$. 
Lemma. Let $\mu$ be any finite measure on $\mathbb R$ with finite first $k+1$ moments and with support set $S_\mu$ of cardinality at least $k+2$, where $k$ is a natural number. For any interval $J\subseteq\mathbb R$, let 
$N_k(\mu,J)$ denote the set of all measures $\nu$ on $\mathbb R$ such that $\int_J x^j\nu(dx)=\int_J x^j\mu(dx)$ for all $j=0,\dots,k$, with $N_k(\mu):=N_k(\mu,\mathbb R)$. Let $M_k(\mu):=\sup\{\int_{\mathbb R}x^{k+1}\nu(dx)\colon\nu\in N_k(\mu)\}$ and 
$m_k(\mu):=\inf\{\int_{\mathbb R}x^{k+1}\nu(dx)\colon\nu\in N_k(\mu)\}$. Then 
$$m_k(\mu)<\int_{\mathbb R}x^{k+1}\mu(dx)<M_k(\mu),\tag{-1+}
$$ 
with an emphasis on the strictness of the inequalities.
Proof. To obtain a contradiction, suppose that $\int_{\mathbb R}x^{k+1}\mu(dx)=M_k(\mu)$. 
Since $S_\mu$ is of cardinality at least $k+2$, there are pairwise disjoint closed intervals $J_0,\dots,J_{k+1}$ of nonzero $\mu$-measure. Again by 
my answer to the question at 
[mathoverflow.net/questions/229646], 
for each $i$ there is a measure $\nu_i\in N_{k+1}(\mu,J_i)$ with a finite support set contained in $J_i$, so that $p_i:=\nu_i(\{x_i\})>0$ for some $x_i\in J_i$. It also follows that the points $x_i$ are pairwise distinct. 
Let now $\nu(A):=\mu(A\setminus J_0\setminus\dots\setminus J_{k+1})+\sum_{i=0}^{k+1}\nu_i(A)$ for all Borel sets $A\subseteq\mathbb R$. 
Then $\nu\in N_{k+1}(\mu)$ and, in particular, $\int_{\mathbb R}x^{k+1}\nu(dx)=M_k(\mu)$. 
For the pairwise distinct points $x_i$ mentioned above, consider the system of equations 
$$\sum_{i=0}^{k+1}a_i x_i^j=I\{j=k+1\}\quad\text{for } j=0,\dots,k+1,
$$
for the unknowns $a_0,\dots,a_{k+1}$, where $I\{\cdot\}$ is the indicator function. The determinant of this linear system is a nonzero Vandermonde determinant, and so, the system has a solution $(a_0,\dots,a_{k+1})$. For that solution and real $t>0$, let $\nu_t:=\nu+t\sum_{i=0}^{k+1}a_i\delta_{x_i}$, where $\delta_x$ is the Dirac measure at $x$. 
If $t$ is small enough (so that $p_i+ta_i\ge0$ for all $i$), then $\nu_t\in N_k(\mu)$ and $\int_{\mathbb R}x^{k+1}\nu_t(dx)=\int_{\mathbb R}x^{k+1}\nu(dx)+t=M_k(\mu)+t>M_k(\mu)$, which is a contradiction. 
Similarly, the assumption $\int_{\mathbb R}x^{k+1}\mu(dx)=m_k(\mu)$ leads to a contradiction. QED
A: If by "a Gaussian mixture distribution" you mean the mixture of zero-mean Gaussian distributions (as your example may suggest), then the answer is no, for any $\epsilon\ne0$ and any $k\ge2$ -- except for the trivial case $k=2$ and the case $k=4$ addressed in your example. Indeed, by the symmetry, without loss of generality $k$ is an even number $\ge6$, so that $m:=k/2$ is a natural number $\ge3$. The random variable (r.v.) $X$ is equal to $SZ$ in distribution, where $S$ is a nonnegative r.v. independent of $Z$. At that, $E(SZ)^2=EZ^2$ and $E(SZ)^4=EZ^4$, whence $ES^2=ES^4=1$, which implies that $S^2=1$ (and hence $S=1$) almost surely, so that $EX^j=E(SZ)^j=EZ^j$ for all natural $j$.
