From the get-go we should be able to tell the equation (2) is incorrect, since it equals
$$\left(\int_{-\infty}^\infty p(x_1)dx_1\right)\cdots\left(\int_{-\infty}^\infty p(x_{n-2})dx_{n-2}\right)\int_{-\infty}^\infty p(x_{n-1})p(x-x_{n-1})dx_{n-1} \tag{2-bad}$$
which is $1\cdots 1\cdot P_2(x)$ instead of the $P_n(x)$ it's supposed to be.
The correct form is something that can be intuited. Since step sizes are independent events, the probability density of having step sizes $(x_1,\cdots,x_n)$ should be $p(x_1)\cdots p(x_{n-1})p(x_n)$ which we should then integrate over the hyperplane defined by the equation $x_1+\cdots+x_n=x$ in order to obtain the probability the $n$ steps sum to $x$. Such a domain of integration can be parametrized by letting the variables $x_1,\cdots,x_{n-1}$ roam freely and then setting $x_n=x-(x_1+\cdots+x_{n-1})$, yielding
$$\int_{-\infty}^\infty\cdots\int_{-\infty}^\infty p(x_1)\cdots p(x_{n-1})p\big(x-(x_1+\cdots+x_{n-1})\big) dx_{n-1}\cdots x_1. \tag{2-good}$$
We can derive this from (1) using substitution if we wish. Observe
$$P_n(x)=\int_{-\infty}^\infty p(x-u_1)P_{n-1}(u_1)du_1=\int_{-\infty}^\infty\int_{-\infty}^\infty p(x-u_1)p(u_1-u_2)P_{n-2}(u_2)du_2du_1$$
$$\cdots =\int_{-\infty}^\infty\cdots\int_{-\infty}^\infty p(x-u_1)p(u_1-u_2)\cdots p(u_{n-2}-u_{n-1})P_1(u_{n-1})du_{n-1}\cdots du_1 $$
and of course $P_1=p$. Make a change of variables
$$\begin{bmatrix}x_1 \\ x_2 \\ \vdots \\ x_{n-2} \\ x_{n-1}\end{bmatrix}=\begin{bmatrix}u_1-u_2 \\ u_2-u_3 \\ \vdots \\ u_{n-2}-u_{n-1} \\ u_{n-1}\end{bmatrix}=\begin{bmatrix} 1 & -1 & 0 & \cdots & 0 & 0 \\ 0 & 1 & -1 & \cdots & 0 & 0 \\ 0 & 0 & 1 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 1 & -1 \\ 0 & 0 & 0 & \cdots & 0 & 1 \end{bmatrix} \begin{bmatrix}u_1 \\ u_2 \\ \vdots \\ u_{n-2} \\ u_{n-1}\end{bmatrix}.$$
The Jacobian determinant of this upper triangular matrix is $1$, and $u_1=x_1+\cdots+x_{n-1}$, so the integral becomes precisely the one written in $(\text{2-good})$.
I'm not aware of a good reference for the formula I gave; here's an explanation instead, (albeit a somewhat hand-waving one).
If integrals like the one you gave are well defined, multiplying by $\delta(x_i)$ where $x_i$ is a coordinate ought be the same as setting those coordinates to zero:
$$
\int_{U\subseteq\mathbb{R}^3}f(x_1,x_2,x_3)\delta(x_2)\delta(x_3)dx_1dx_2dx_3=\int_{\{x_1\in\mathbb{R}:(x_1,0,0)\in U\}}f(x_1,0,0)dx_1
$$
and so on for other dimensions and other numbers of deltas. Additionally, we would expect such integrals to still obey the change of variables formula:
$$
\int_Vf(\mathbf{x})d\mathbf{x}=\int_Uf(\varphi(\mathbf{y}))|\det(D\varphi(\mathbf{y}))|d\mathbf{y}
$$
Where $\varphi:U\to V$ is a diffeomorphism between open subsets of $\mathbb{R}^n$.
If we assume that the zero sets of $g$ and $h$ are transverse (i.e. $\nabla g$ and $\nabla h$ are linearly independent on $g^{-1}(0)\cap h^{-1}(0)$), then this is already enough. Let $\mathbf{p}\in g^{-1}(0)\cap h^{-1}(0)$. By inverse function theorem there exist a set of local coordinates $\mathbf{y}:V\to U$ on a neighborhood $V$ of $\mathbf{p}$ such that $\mathbf{y}(p)=\mathbf{0}$, $y_3=g$, and $y_4=h$. Let $\varphi:U\to V$ be the inverse of these coordinates. Applying the change of variables formula gives
$$
\int_Vf(\mathbf{x})\delta(g(\mathbf{x}))\delta(h(\mathbf{x}))d\mathbf{x}=\int_Uf(\varphi(\mathbf{y}))\delta(g(\varphi(\mathbf{y}
))\delta(h(\varphi(\mathbf{y}))|\det(D\varphi(\mathbf{y})|d\mathbf{y}
$$
Since $\varphi$ is inverse to $\mathbf{y}$, we can rewrite this as
$$
=\int_U\frac{f(\varphi(\mathbf{y}))\delta(y_3)\delta(y_4)}{|\det(D\mathbf{y}(\varphi(\mathbf{y})))|}d\mathbf{y}
$$
Now that the deltas are composed with coordinates, we can restrict to a surface integral.
$$
=\int_{U\cap\mathbb{R}^2}\frac{f(\varphi(y_1,y_2,0,0))}{|\det(D\mathbf{y}(\varphi(y_1,y_2,0,0)))|}dy_1dy_2
$$
The columns of $D\mathbf{y}$ are just the gradients of the coordinates.
$$
=\int_{U\cap\mathbb{R}^2}\frac{f(\varphi(\mathbf{y}))\delta(y_3)\delta(y_4)}{|\det([\nabla y_1,\nabla y_2,\nabla g,\nabla h])|(\varphi(y_1,y_2,0,0))}dy_1dy_2
$$
Let $\pi$ be the orthogonal projection onto $T(g^{-1}(0)\cap h^{-1}(0))$. Since this amounts to subtracting multiples of $\nabla g$ and $\nabla h$, we can apply it to $\nabla y_1$ and $\nabla y_2$ without changing the determinant:
$$
=\int_{U\cap\mathbb{R}^2}\frac{f(\varphi(y_1,y_2,0,0))}{|\det([\pi(\nabla y_1),\pi(\nabla y_2),\nabla g,\nabla h])|(\varphi(y_1,y_2,0,0))}dy_1dy_2
$$
To split the determinant, note that $|\det([a,b,c,d])|$ is equal to the volume of the parallelepiped spanned by $a,b,c,d$. If $\operatorname{span}(a,b)\perp\operatorname{span}(c,d)$, then this is equal to the product of areas the parallelograms spanned by $(a,b)$ and $(c,d)$. In terms of wedge products, this gives
$$
=\int_{U\cap\mathbb{R}^2}\frac{f(\varphi(y_1,y_2,0,0))}{\|\nabla g\wedge\nabla h\|(\varphi(y_1,y_2,0,0))}\frac{dy_1dy_2}{\|\pi(\nabla y_1)\wedge\pi(\nabla y_2)\|(\varphi(y_1,y_2,0,0))}
$$
Where the last term is just the standard area element of $g^{-1}(0)\cap h^{-1}(0)$ in the coordinates $y_1,y_2$. Writing this in a more coordinate-independent way, we (almost) have the desired formula.
$$
=\int_{V\cap(g^{-1}(0)\cap h^{-1}(0))}\frac{f}{\|\nabla g\wedge\nabla h\|}dA
$$
This result is only local, but you can recover the global version using a partition of unity.
Best Answer
You can compute the Fourier transform of $f$: $$\widehat{f}(\xi)=\int f(x) e^{ix \xi}dx=\int \int\cdots\int\delta\left(\sum_{n=1}^{N}x_{n}-x\right)\prod_{n=1}^{N}\frac{1}{1+x_{n}^2} e^{i x \xi}dx_1\cdots dx_N dx$$ so: \begin{align}\widehat{f}(\xi)&= \int\cdots \int\left( \int\delta\left(\sum_{n=1}^{N}x_{n}-x\right) e^{i x \xi} dx \right)\prod_{n=1}^{N}\frac{1}{1+x_{n}^2} dx_1\cdots dx_N\\ &= \int\cdots \int\left( \prod_{n=1}^{N} e^{i x_n \xi} \right)\prod_{n=1}^{N}\frac{1}{1+x_{n}^2} dx_1\cdots dx_N\\ &=\left( \int e^{i x \xi} \frac{1}{1+x^2}\right)^N=\pi^N e^{-N |x|} \end{align} and using the inverse Fourier transform one can obtain your result.