Recently I was trying to show CLT by analyzing convolution of pdfs of infinitely large number of i.i.d random variables. I'm saying "show" instead of "prove" because I don't have enough mathematical experience to produce rigorous proofs. A case with CLT was ok, but when I tried to use the same method to get something like Berry-Esseen theorem I stumbled upon a problem that I couldn’t solve.
So, suppose $X_1, X_2, ..., X_n$ are i.i.d random variables with $\mathbb{E}X=0$ and $ \mathbb{Var}X=1$.
We want to find limit of distribution $p_{S_n}(t)$ of $S_n=\frac{X_1+X_2+...+X_n}{\sqrt{n}}$ as $n\to\infty$. After adding one more $X_i$ we will have a new random variable $S_{n+1}$ such that
$$S_{n+1} = \frac{\sqrt{n}}{\sqrt{n+1}}S_n+\frac{X_{n+1}}{\sqrt{n}}$$
Denoting $\frac{X_{n+1}}{\sqrt{n}}=\xi$ with pdf $p_{\xi}(x)$, we have
$$p_{S_{n+1}}(t) = \int p_{\xi}(x) \frac{\sqrt{n+1}}{\sqrt{n}}p_{S_n}\Big(\frac{\sqrt{n+1}}{\sqrt{n}}(t-x)\Big)dx$$
Some simplifications can be made. First, $n \gg 1$, so $\frac{\sqrt{n+1}}{\sqrt{n}} \approx 1+\frac{1}{2n}$. Second, $\xi$ has very small variance of $\frac{1}{n+1}$, so all its mass is in the vicinity of 0 . So, integral is nonzero only for sufficiently small values of x. And that's why we can expand $p_{S_n}$. $$p_{S_n}\Big(\frac{\sqrt{n+1}}{\sqrt{n}}(t-x)\Big) \approx p_{S_n}\Big((1+\frac{1}{2n})t\Big)-xp'_{S_n}\Big((1+\frac{1}{2n})t\Big)+\frac{x^2}{2}p''_{S_n}\Big((1+\frac{1}{2n})t\Big)-\frac{x^3}{6}p'''_{S_n}\Big((1+\frac{1}{2n})t\Big)+...$$
$\mathbb{E}\xi=0$ and $ \mathbb{Var}\xi\approx\frac{1}{n}$ so first derivative term will be zero after integration and second derivative term will of order $\frac{1}{n}$. And there is no need to expand inside $p''_{S_n}$ and $p'''_{S_n}$ if we want to consider only terms of order bigger than $\frac{1}{n^2}$. So after integration we have $$p_{S_{n+1}}(t)=\Big(1+\frac{1}{2n}\Big)\Big(p_{S_n}\Big((1+\frac{1}{2n})t\Big)+0+\frac{1}{2(n+1)}p''_{S_n}(t)-\frac{\mu_3}{6(n+1)^{\frac{3}{2}}}p'''_{S_n}(t)\Big)\approx\\ \approx p_{S_n}(t)+\frac{p_{S_n}(t)}{2n}+\frac{tp'_{S_n}(t)}{2n}+\frac{p''_{S_n}(t)}{2n}-\frac{\mu_3p'''_{S_n}(t)}{6n^{\frac{3}{2}}} \tag{1} $$
As $n\to\infty$, $p_{S_{n+1}}-p_{S_n}\to0$, so after multiplying by $2n$ we have a differential equation: $$p(t)+tp'(t)+p''(t)-\frac{\mu_3p'''(t)}{3\sqrt{n}}=0$$ If we ignore the last term and take into account that $p(-\infty)=0$, we'll have:
$$(tp(t)+p'(t))'=0$$ $$tp(t)+p'(t)=0$$ $$p(t)=ce^{-\frac{t^2}{2}}$$
And this is essentially CLT. So far so good.
But if we denote $\varepsilon=\frac{\mu_3}{3\sqrt{n}}$ and consider DE with the fourth term included we'll have:
$$p(t)+tp'(t)+p''(t)-\varepsilon p'''(t)=0$$ $$tp(t)+p'(t)-\varepsilon p''(t)=0$$ Let's look for a solution of the form $$p=p_0+\varepsilon p_1+...$$ We have then for different orders of $\varepsilon$: $$\varepsilon^0:tp_0+p'_0=0\Rightarrow p_0=c_0e^{-\frac{t^2}{2}}$$ $$\varepsilon^1:tp_1+p'_1=p''_0\Rightarrow p_1=c_1e^{-\frac{t^2}{2}}+c_0e^{-\frac{t^2}{2}}\big(\frac{t^3}{3}-t\big)$$ So final solution is: $$p(t)=(c_0+\varepsilon c_1)e^{-\frac{t^2}{2}}+\varepsilon c_0 e^{-\frac{t^2}{2}}\big(\frac{t^3}{3}-t\big)$$ And here is the problem. If I determine $c_0$ from the equation of order $\varepsilon^0$ (from $\int p dt=1$), then $c_0=\frac{1}{\sqrt{2\pi}}$ and $c_1=0$. But then on the other hand $\int t^3pdt=\frac{\mu_3}{\sqrt{n}}=3\varepsilon$ and so $$3\varepsilon=\int t^3pdt=\int\frac{t^3}{\sqrt{2\pi}}e^{-\frac{t^2}{2}}+\frac{\varepsilon t^3}{\sqrt{2\pi}} e^{-\frac{t^2}{2}}\big(\frac{t^3}{3}-t\big)dt=0+\frac{\varepsilon}{\sqrt{2\pi}}2\sqrt{2\pi}=2\varepsilon\neq3\varepsilon$$ And if I try to determine $c_0$ and $c_1$ from the full solution, I get $$c_0=\frac{3}{2}, c_1=\frac{1}{\varepsilon}\big(\frac{1}{\sqrt{2\pi}}-\frac{3}{2}\big)$$ The fact that $c_1\to\infty$ as $\varepsilon\to0$ doesn't seem right for many reasons.
So as far as I understand, the cause of these contradictions is that multiplier of $\frac{1}{\sqrt{n}}$ that I got is $\frac{\mu_3}{3}$ instead of $\frac{\mu_3}{2}$. I have my suspicions about where things might have gone wrong, but I could not find a convincing way to fix my mistake.
For example, one way to get the "right" answer is to modify equation (1) and consider $$p_{S_{n+2}}-p_{S_{n+1}}\approx p_{S_{n+2}} - p_{S_n}-\frac{p_{S_n}}{2n}-\frac{tp'_{S_n}}{2n}-\frac{p''_{S_n}}{2n}+\frac{\mu_3p'''_{S_n}}{6n^{\frac{3}{2}}}\approx\\ \approx p_{S_{n+1}} + \frac{p_{S_{n+1}}}{2(n+1)}+\frac{tp'_{S_{n+1}}}{2(n+1)}+\frac{p''_{S_{n+1}}}{2(n+1)}-\frac{\mu_3p'''_{S_{n+1}}}{6(n+1)^{\frac{3}{2}}} - p_{S_n}-\frac{p_{S_n}}{2n}-\frac{tp'_{S_n}}{2n}-\frac{p''_{S_n}}{2n}+\frac{\mu_3p'''_{S_n}}{6n^{\frac{3}{2}}}\approx\\ \approx p_{S_{n+1}} - p_{S_{n}} -\frac{p_{S_n}}{2n^2}-\frac{tp'_{S_n}}{2n^2}-\frac{p''_{S_n}}{2n^2}+\frac{\mu_3p'''_{S_n}}{4n^{\frac{5}{2}}}$$ And so as $(p_{S_{n+2}}-p_{S_{n+1}}) - (p_{S_{n+1}}-p_{S_{n}})\to0$ we get an equation $$\frac{\mu_3}{2\sqrt{n}} p''' - p''-p'-tp=0$$ Its solution (assuming p is continuous) with $c_0=\frac{1}{\sqrt{2\pi}}$ after integrating from $-\infty$ to x corresponds to result of Berry and Esseen namely $$F_n(x) \approx \Phi(x) + \frac{1}{\sqrt{n}}\frac{\mu_3}{6\sqrt{2\pi}}e^{-\frac{x^2}{2}}\big(1-x^2\big)$$ where $F_n$ is a CDF of $p_{S_n}$ and $\Phi$ is a CDF of standard Normal.
But I don't know why would I want to consider second differences in the first place. And moreover, it seems that for instance, third differences give wrong result again.
So what is a proper way to get a right result with this approach? I would appreciate if you could help me with that.