Moment matching: construction of a mixture of Gaussian distribution with lower moments identical to Gaussian

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