Moment method for weak convergence 
Say, we have a sequence of probability distributions $\mu_n$ on $\mathbb R$, that are uniformly subgaussian in the sense that $$\mu_n(\mathbb R\setminus[-R,R])\leq Ce^{-CR^2}$$ for some positive constant $C$. Assume that all moments of these distributions exist and that they converge to those of a standard normal distribution, i.e., $$\lim_{n\to\infty}\int x^k~d\mu_n(x)=\begin{cases}(k-1)!!&\text{if }k\text{ is even}\\ 0 &\text{else}\end{cases}.$$ Then $\mu_n$ converges weakly to a standard-normal distribution $\mu$, i.e., $$\lim_{n\to\infty} \int f~d\mu_n=\int f~d\mu$$ for all bounded continuous $f$. 

The above claim is clearly true: The uniform subgaussian hypothesis implies that the $k$-th moments are uniformly bounded by $(Ck)^{k/2}$ and therefore the characteristic functions can be written as $$\int e^{itx}~d\mu_n(x)=\sum_{k=0}^\infty \frac{(it)^k}{k!}\int x^k~d\mu_n(x),$$ which converges point wise to the characteristic function of the normal distribution. The result now follows from Levy's Continuity Theorem.
Question: Is there a way to prove the above statement directly? (I am concerned with something related where the measures are random themselves and I can't argue with the characteristic function)
Thoughts: Fix $\epsilon>0$ and $f\in C_b(\mathbb R)$ and choose $\lambda<\infty$ such that $\mu_n(\mathbb R\setminus[-\lambda,\lambda])\leq \epsilon$. From the Weierstraß approximation theorem, we find a polynomial $p$ that approximates $f$ on the interval $[-\lambda,\lambda]$ up to an error of $\epsilon$. Thus we find 
$$\big\lvert\int f~d\mu_n-\int f~d\mu\big\lvert\leq \int_{[-\lambda,\lambda]}\lvert f-p\lvert~d\mu_n+\big\lvert\int_{[-\lambda,\lambda]^c} f-p~d\mu_n\big\lvert+\big\lvert\int p~d(\mu_n-\mu)\big\lvert+\int \lvert p-f\lvert~d\mu.$$ The first term is bounded by $\epsilon$ and the third term converges to zero by assumption. The fourth term can be bound by $\epsilon$ and the integral outside the finite interval. Since $f$ is bounded we can estimate $\lvert f(x)-p(x)\lvert\leq c x^{2k}$ outside the interval for some $c,k\in\mathbb N$. It would therefore remain to show that $$\int_{[-\lambda,\lambda]^c} x^{2k}~d(\mu_n+\mu)$$ is small. Does this somehow follow from the subgaussian assumption (which wasn't really used till now)? 
 A: In the proof of the claim by characteristic function, we have to check that we can apply the dominated convergence theorem, which can be done using the assumption of uniform subgaussinity. 
Note that the measure $\mu$ is Gaussian, hence it has moments of any order, which implies 
$$\int_{\mathbf R\setminus [-\lambda,\lambda]}  x^{2 j}\mathrm d\mu \leqslant \frac 1{\lambda^2}\int_{\mathbf R}  x^{2j+2}\mathrm d\mu $$
and this can be made arbitrarily small. 
Notice that if $\lambda$ is an integer, then 
$$\int_{\mathbf R\setminus [-\lambda,\lambda]}  x^{2 j}\mathrm d\mu_n 
=\sum_{l\geqslant\lambda}\int_{|x|\in[l,l +1)}x^{2 j}\mathrm d\mu_n\leqslant \sum_{l\geqslant\lambda}l^{2j} \int_{|x|\geqslant l}\mathrm d\mu_n  . $$
Using the assumption of subgaussianity, we infer the bound 
$$\int_{\mathbf R\setminus [-\lambda,\lambda]}  x^{2 j}\mathrm d\mu_n\leqslant C\sum_{l\geqslant\lambda}l^{2j}e^{-Cl ^2}.$$
But the problem is the following: the chosen $\lambda$ depends on $j$, which depends itself of $p$, which depends... on $\lambda$. 
Let $f$ be a continuous and bounded function. Define for a fixed $R$ and and integer $l$ the polynomial 
$$P_{l,R}(x)=\sum_{j=0}^l\binom jl f\left(\frac{2R}lj-R\right)
\left(\frac{x+R}{2R}\right)^j
\left(\frac{R-x}{2R}\right)^{l-j}.$$
Then for each $\delta$ and each $x\in [-R,R]$, the following estimate holds 
$$|f(x)-P_{l,R}(x)|\leqslant \frac 2{l\delta^2}\sup_{t\in [-R,R]}|f(t)|R^2
+\sup_{t,s\in [-R,R],|s-t|\lt\delta}|f(t)-f(s)|.$$
This follows from the estimate we obtain when we consider the unit interval instead of $[-R,R]$, which are obtained by splitting the sum defining $P_{l,R}$ over the $j$ for which $|j/l-x|$ is smaller of greater than $\delta$.
We thus have 
\begin{align}
\left|\int fd\mu_n-\int fd\mu\right|&\leqslant \frac 2{l\delta^2}\sup_{t\in [-R,R]}|f(t)|R^2
+\sup_{t,s\in [-R,R],|s-t|\lt\delta}|f(t)-f(s)|\\
&+\lVert f\rVert_\infty\left(\mu_n(\mathbf R\setminus [-R,R]\right)+\mu\left(\mathbf R\setminus [-R,R]\right)+\\
&+\left|\int P_{l,R}d\mu_n-\int P_{l,R}d\mu\right|+\left|\int_{\mathbf R\setminus [-R,R]} P_{l,R}d\mu_n-\int_{\mathbf R\setminus [-R,R]} P_{l,R}d\mu\right|.
\end{align}
Using the assumption of convergence of the moments, we infer the upper bound 
\begin{align}
\limsup_{n\to\infty}\left|\int fd\mu_n-\int fd\mu\right|&\leqslant \frac 2{l\delta^2}\sup_{t\in [-R,R]}|f(t)|R^2
+\sup_{t,s\in [-R,R],|s-t|\lt\delta}|f(t)-f(s)|\\
&+\lVert f\rVert_\infty\left(\sup_n\mu_n(\mathbf R\setminus [-R,R]\right)+\mu\left(\mathbf R\setminus [-R,R]\right)+\\
&+\limsup_{n\to\infty}\left|\int_{\mathbf R\setminus [-R,R]} P_{l,R}d\mu_n\right|+\left|\int_{\mathbf R\setminus [-R,R]} P_{l,R}d\mu\right|.
\end{align}
Notice that 
$$|P_{l,R}(x)|\leqslant \lVert f\rVert_\infty\left(\frac{|x+R|+|x-R|}{2R}\right)^l,$$
hence 
$$\int_{[R,\infty)}|P_{l,R}(x)| d\mu_nx\leqslant \lVert f\rVert_\infty R^{-l}
\int_{[R,\infty)}x^ld\mu_nx\leqslant C\sum_{j\geqslant R}(j+1)^lCe^{-Cj^2}R^{-l},$$
which can be bounded by 
$$l!e^{R-CR^2}\sum_j e^{j-Cj^2}.$$
We can thus choose $l:=R^3$.
