Is a probability density function necessarily a $L^2$ function?

If a nonnegative continuous real valued function $f$ is integrable over $\mathbb{R}$ with $$\int_\mathbb{R} f\,\mathrm{d}x = 1,$$ does it hold true $$\int_\mathbb{R} f^2 \,\mathrm{d}x<\infty?$$

Motivation: I am wondering if the mean squared error (MSE) is well defined, since we need the target density function to be in $L^2$ space in order to have a finite MSE.

The answer is no. But your continuity requirement makes it more difficult to find an explicit (constructed) example.

Define the function $\varphi$ on $\mathbb{R}$ as: $$\forall x\in\mathbb{R},\quad\varphi(x)=\begin{cases}\sin(x)&\text{if $x\in[0,\pi]$}\\0&\text{otherwise.}\end{cases}$$ Clearly, the function $\varphi$ is continuous on $\mathbb{R}$ (it's just a single arch of sine). It is well-known that $\lVert\varphi\rVert_1=2$ and that $\lVert\varphi\rVert_2^2=\pi/2$.

Define the sequence of function $(\varphi_n)_{n\geq1}$ on $\mathbb{R}$ by: $$\forall n\geq1,\quad\forall x\in\mathbb{R},\quad\varphi_n(x)=\varphi\bigl(8^n(x-2n\pi)\bigr).$$

It should be clear that the support of the $\varphi_n$'s are all distinct: the support of $\varphi_n$ is $(2n\pi,8^{-n}\pi+2n\pi)$.

Moreover, for all $n\geq1$, $$ \lVert\varphi_n\rVert_1=\int_{\mathbb{R}}\varphi_n(x)\,\mathrm{d}x=\frac2{8^n}\quad\text{and}\quad\lVert\varphi_n\rVert_2^2=\int_{\mathbb{R}}\varphi_n(x)^2\,\mathrm{d}x=8^{-n}\frac{\pi}2.$$

Now, define the series of functions $f$ by: $$\forall x\in\mathbb{R},\quad f(x)=\frac12\sum_{n=1}^{+\infty}4^n\varphi_n(x)$$ (the graph of $f$ consists of bumps that are narrower but higher). It should be clear that $f$ is well-defined and continuous and non-negative. Since the support of the $\varphi_n$'s are distinct, we can integrate $f$ term by term: $$\lVert f\rVert_1=\int_{\mathbb{R}}f(x)\,\mathrm{d}x=\frac12\sum_{n=1}^{+\infty}4^n\int_{\mathbb{R}}\varphi_n(x)\,\mathrm{d}x=\sum_{n=1}^{+\infty}\frac{4^n}{8^n}=1.$$

Also, since the support of the $\varphi_n$'s are distinct, for all $x\in\mathbb{R}$: $$f(x)^2=\frac14\sum_{n=1}^{+\infty}16^n\varphi_n(x)^2$$ and $$\lVert f\rVert_2^2=\int_{\mathbb{R}}f(x)^2\,\mathrm{d}x =\frac14\sum_{n=1}^{+\infty}16^n\lVert\varphi_n\rVert_2^2=\frac14\sum_{n=1}^{+\infty}2^n\frac\pi2=+\infty.$$