The "correct" way to approximate $\text{var}(f(X))$ via Taylor expansion tl;dr: There are two commonly reported formulas for approximating $\text{var}(f(X))$, but one is notably better than the other. Since it isn't the "standard" Taylor expansion, where does it come from, and why is it better?
Details: Let $X$ be a real random variable and $f:\mathbb{R}\to\mathbb{R}$. There is a standard way to approximate $\text{var}(f(X))$ using Taylor expansions of moments, e.g. $E f(X)$, $E f^2(X)$, etc. Doing so yields the following second-order approximation:
\begin{align}
\operatorname{var}(f(X))\approx [f'(EX)]^2\operatorname{var}(X)-\frac{[f''(EX)]^2}{4}\operatorname{var}^2(X)
:= V_1.
\end{align}
For a formal proof, see this post and answer. Oddly, this is not the formula quoted in the corresponding Wikipedia page:
\begin{align}
\operatorname{var}(f(X))\approx [f'(EX)]^2\operatorname{var}(X) + \frac{[f''(EX)]^2}{2}\operatorname{var}^2(X)
:= V_2.
\end{align}
Note the difference in the coefficient of second term: $-1/4$ vs $+1/2$.
At first, I assumed this was a typo in the Wikipedia page. After running some quick simulations, however, it seems that the second approximation $V_2$ is much better than $V_1$! (Admittedly I did not run exhaustive tests but after a few dozen examples the difference was quite clear.) 
In fact, here is a partial explanation for why $V_1$ often fails catastrophically: If $f'(EX)\approx 0$, then $V_1 < 0$. Presumably, this can be corrected by using a third-order approximation.
My question, though, is (a) How do we derive $V_2$, and (b) Why does it outperform $V_1$?
 A: I cannot speak to the derivation of the first approximation (which looks wrong to me).  However, the second equation is obtained using a second-order Taylor approximation to $f$ for the case where the underlying distribution is centred, unskewed and mesokurtic.  In this case, you have $\mu=0$, $\gamma=0$ and $\kappa=3$.  Using the general form of the Taylor approximation you obtain:
$$\begin{equation} \begin{aligned}
\mathbb{V}[f(X)]
&\approx ( f''(\mu)^2 \mu^2 - f'(\mu)f''(\mu) \mu + f'(\mu)^2 ) \cdot \sigma^2 \\[6pt]
&\quad - \frac{f''(\mu)(f'(\mu) + \mu f''(\mu))}{2} \cdot \gamma \sigma^3 
+ \frac{f''(\mu)^2}{4} \cdot (\kappa-1) \sigma^4 \\[6pt]
&= f'(\mu)^2 \cdot \sigma^2 + \frac{f''(\mu)^2}{2} \cdot \sigma^4. \\[6pt]
\end{aligned} \end{equation}$$
The first approximation does not look correct to me, and I see no evidence that it is a "commonly reported formula".  This approximation cannot be derived from the general second-order Taylor approximation for any assumed level of kurtosis, so I find it unsurprising that it performs poorly.  (It would require $\kappa = 0$ which is not a valid kurtosis value.)  For this reason, I would expect the second approximation to perform better than the first, except possibly in the case where the kurtosis of the underlying distribution is highly platykurtic.
