You've correctly identified your problem: The fact that $f$ is analytic doesn't imply that $f(z) = z$ (certainly $e^z$, $z^2$, $\cos z$ are analytic). Of course, if you were coming up with a counterexample, it would be perfectly fine to let $f(z) = z$ (or whatever the counterexample should be) and prove the statement false. However, the statement actually is true, so this approach is hopeless.
One way to do this is to use the fact that $f$ is analytic if it's equal to its Taylor series. Now using the continuity of the operation of complex conjugation is analytic (together with the fact that $\overline z^k = \overline{z^k}$), so if
$$f(z) = \sum_{k = 0}^{\infty} a_k z^k$$
we have
\begin{align*}
\overline{f({\overline z})} &= \overline{\sum_{k = 0}^{\infty} a_k \overline z^k} \\
&= \sum_{k = 0}^{\infty} \overline{a_k \overline z^k} \\
&= \sum_{k = 0}^{\infty} \overline{a_k} z^k
\end{align*}
Now since $|a_k| = |\overline{a_k}|$, both series have the same (infinite) radius of convergence, so $\overline{f(\overline z)}$ is simultaneously analytic with $f$. Now the square of an analytic function is analytic, and we're done.
Okay, on closer sketching, this isn't actually that bad. For the sake of convenience, let's assume instead that $Z\sim \mathcal{N}(0,I)$ when viewed as an $\mathbb{R}^2$-valued random variable. Denote its density $f$.
Define $P:\mathbb{R}^2\setminus \{x\in (0,\infty),y=0\}\to (0,\infty)\times (0,2\pi)$ to be the standard polar coordinate transformation. Then, since $\{x\in (0,\infty),y=0\}$ is a $\mathcal{N}(0,I)$-null set, we can apply the Jacobi Coordinate Transformation theorem to get that $(R,\Theta):=P(Z)$ has density
$$rf(P^{-1}(r,\theta))=rf(r(\cos(\theta)+i\sin(\theta))=\frac{r}{2\pi}\exp(-\frac{r^2}{2})=\frac{1}{2\pi}\cdot r\exp(-\frac{r^2}{2}),$$
which is clearly a factorisation of the density, implying that $R$ and $\Theta$ are independent, and $\Theta$ is uniformly distributed on $(0,2\pi)$. Note that $R$ and $\Theta$ clearly have moments of all orders.
Accordingly, we get, by applying independence coordinate-wise, that
$$
E Z^m \overline{Z^k}=E(R^{m+k} e^{i (m-k)\Theta})=E(R^{m+k})E(e^{i(m-k)\Theta}),
$$
and
$$
E(e^{i(m-k)\Theta})=\frac{1}{2\pi}\left(\int_0^{2\pi} \cos((m-k)\theta)\textrm{d}\theta+i\int_0^{2\pi}\sin((m-k)\theta)\textrm{d}\theta\right)=0,
$$
since $m\neq k$. This yields the desired.
Best Answer
When you're reviewing for an entire course, you probably know that a function cannot have an antiderivative in a region unless it is differentiable there. So it suffices to show that your $f$ is not differentiable in any nonempty open set.
A straightforward approach would be to show that it doesn't satisfy the Cauchy-Riemann equations anywhere except at the origin -- that's not a lot of computation in this case.
Alternatively you could observe that $\bar z^2 = \overline{z^2}$, if you have proved a theorem that the only way for a function and its complex conjugate both to be differentiable is if it is constant, which $z^2$ isn't. (You may know this theorem in the form: If $g$ is real-valued and differentiable everywhere in a connected open set, then it is constant. This can be applied to $g(z)=f(z)+\overline{f(z)}$ and then $g_2(z)=\frac{f(z)-\overline{f(z)}}i$).