Variance and asymptotic normality of $\frac{1}{n-1}\sum_{i=1}^{n-1}(x_{i+1}-x_i)^2$, where $X \sim \mathcal{N}(0,1)$ Consider a length $n$ vector $\mathbf{x}$ containing $n$ i.i.d. observations $\{x_i\}_{i=1}^n$ of a standard normal random variable $X$. Let $\mathbf{z}$ be a length $n-1$ vector whose entries are $z_i = x_{i+1}-x_i$. I have a statistic $s(\mathbf{z}) = \frac{1}{n-1}\sum_{i=1}^{n-1}(z_i)^2$ that is meant to measure how close to sorted a given permutation of random samples is, and I want to know its properties under the null hypothesis where the order is random, for the purpose of calculating a $p$-value.
I have two questions:


*

*Is $s(\mathbf{z})$ asymptotically normal?

*What is the variance of $s(\mathbf{z})$  ?


Both of these questions would be easy if the entries of $\mathbf{z}$ were i.i.d., but I have been getting stuck with the fact that they are dependent. For instance, the value of $x_{i+1}-x_i$ depends on the value of $x_{j+1}-x_j$ since $i \neq j$.
For 1., the answer seems to be "yes" based on simulations, but I'm having trouble proving this.
For 2., I thought of estimating the variance by calculating the sample variance of the squared pairwise Euclidean distances of the entries of $\mathbf{x}$, which would work if each $z_i$ were a random sample drawn from $\{x_i - x_j\}_{i,j \in \{1,\ldots,n\}}$, but this is not the case due to the dependence mentioned above. I have found empirically from simulations that dividing this pairwise distance-based estimate by 2 gives the correct result within simulation error, but this might be a lucky coincidence.
 A: TLDR; $s(z)$ is asymptotically normal, and its variance is $\frac {12} {n-1}$ according to CLT for Markov chains. It can be shown that the distribution is a special case of generalized $\chi^2$ distribution.
Markov Chain approach, asymptotics and variance
The sequence $z_i$ is Markov chain because once you know $z_i$ the value of $z_{i+1}$ doesn't depend on $z_k$ where $k<i$. Therefore, Markov chain CLT is applicable. Here's how we apply it.
The sum, or any linear combination of normal r.v.s, is a r.v. itself. Knowing that $z_i\sim\mathcal N(0,2)$ or $z_i\sim\sqrt 2\space \mathcal N(0,1)$, we know that $z_i^2\sim 2\space\chi^2_1$, see the definition of $\chi^2$ distribution. Thus, $\sigma_z^2=\operatorname{var}[z_i^2]=2^2\times 2=8$.
Markov chain CLT states:
$$\sqrt{n-1}(s(z)-\mu)\sim\mathcal N(0,\sigma^2),$$
where $\mu=E[z_i^2]$ and $\sigma^2  = \sigma_z^2 + 2\sum_{k=1}^\infty \operatorname{cov}( z_{1}^2, z_{1+k}^2)=8+2\times 2=12$. 
Hence $\operatorname{var}[s(z)]=\frac{12}{n-1}$
Here's the proof by simulation (Python):
import numpy as np
n = 51
s = np.mean(np.diff(np.random.randn(10000,n))**2,axis=1)
vars = np.var(s)
print(vars)
print(12/(n-1))

Output:
0.23526746023519335
0.24

Note, that if $z_i^2$ weren't correlated then $s(z)$ would have been from scaled $\chi^2$ distribution with variance $\frac 8 {n-1}$. However, due to  overlapping terms $x_i$ in $z_i$ and $z_{i+1}$ we had to apply modified CLT to obtain the asymptotic distribution of $s(z)$.
Acknowledgements: 
My initial answer, which I updated a few times already, did not account for correlation, which was pointed out by @Sextus Empiricus. Also, I used this answer for $\operatorname{cov}( z_{1}^2, z_{1+k}^2)$, where the correlation $\rho=\operatorname{corr}[z_i,z_{i+1}]=-1/2$ and we know that correlation disappears between $z_i$ and $z_j$ when $|i-j|>1$.
The distribution
Let's start with independent row vector of randoms $X'=(x_1,\dots,x_n)$. We get the row vector of differences $Z'$ aby pplying Toeplitz matrix $B'$ as follows $Z'=X'B'$, where
$$B' = \begin{bmatrix} {-1} & 0&\dots & 0 &0\\
 1 & -1 & \dots&0 & 0  \\
 0 & {1}& \dots & 0& 0\\
 \vdots & \vdots & \vdots & \vdots & \vdots \\
 0& 0 & \dots & 1 & -1 \\
 0& 0 & \dots & 0 & {1}
\end{bmatrix}$$
Your quantity then is a quadratic form 
$$s(z)=\frac 1 {n-1} X'B'BX$$
where $B'B$ has a form of a tridiagonal Toeplitz matrix:

Let's apply eigen decomposition $B'B=P'\Lambda P$ then we have:
$$s(z)=\frac 1 {n-1} X'P'\Lambda PX=\frac 1 {n-1} Y'\Lambda Y$$
where $Y=PX\sim\mathcal N(0,I_{n-1})$, i.e. each $Y_i$ (principal component) is an independent normal.
Hence, $$s(z)=\frac 1 {n-1} \sum_{i=1}^{n-1}\lambda_i Y_i^2$$
where $Y_i^2\sim\chi^2_1$ and $\lambda_i$ are eigenvalues. The eigenvalues of Toeplitz tridiagonal matrices are known to form a sine wave and easy to find, see "The Eigenproblem of a Tridiagonal P-Toeplitz Matrix" by Gover .
So the distribution can be seen as a linear combination of $\chi^2$ variables or a generalized $\chi^2$ distribution.
Miscellaneous
We can define a row vector $V'=(x_1,z_1,\dots,z_{n-1})$, then it can be obtained by applying a matrix $D'$ to the original observations $V'=X'D'$, the matrix $B'$ above is the subset of columns of $D'$:

The matrix $D'D$ looks like this:

We can get the matrix $U'$ that recovers the original vector from $V$ as follows: $X'=V'U'$, and $U'=D'^{-1}$.
The matrix $U'$ is upper unit triangular, meaning $u_{ij}=1_{i\ge j}$:

Matrix $A=U'U$, which appears in the quadratic form, has a very interesting form: $a_{ij} = n+1-min(i,j)$, e.g. $n=5$:

A: Small pedantic note: Below I changed the coefficient into $1/\sqrt{n-1}$ otherwise the limiting distribution will be a degenerate distribution (zero variance). In that case one would also need to subtract the mean of the $z_i^2$. That means, only a scaled and shifted sum like $\sum_{i=1}^{n-1} \frac{(x_{i+1}-x_i)^2-1}{\sqrt{n-1}}$ will approach a normal distribution.
Distribution of $s(\mathbf{z})$ as a linear sum of chi-squared variables
The sum $s(\mathbf{z}) = \frac{1}{\sqrt{n-1}}\sum_{i=1}^{n-1}(z_i)^2$ is similarly distributed as the sum $s(\mathbf{y}) = \frac{1}{\sqrt{n-1}}\sum_{i=1}^{n-1}(y_i)^2$ where the $y_i$ are $n-1$ independent normally distributed variables with variance $\lambda_i = 2 + 2 \cos(\frac{i}{n}\pi)$ $$s(\mathbf{z}) \sim  s(\mathbf{y}) \quad \text{where} \quad y_i \sim N\left(0,\lambda_i \right)$$
Consequences:


*

*The variance of $s(\mathbf{z})^2$ is equal to $1/\sqrt{n-1}$ times the sum of the  variances of the individual terms $y_i^2$ (which relate to scaled $\chi_{(1)}^2$ distributions or to gamma distributions). 
For the individual terms which are the squares of normal distributed variables we have 
$$\begin{array}{rcl}
\text{var}(y_i^2) &=& 2 \text{var}(y_i)^2  \\
   &=&2\left( 2 + 2 \cos\left(\frac{i}{n}\pi\right)\right)^2 \\
&=& 12 + 4\cos\left(2\frac{i}{n}\pi\right) + 16 \cos\left(\frac{i}{n}\pi\right) \end{array}$$
and for the sum
$$\begin{array}{rcl}
\text{var}[s(\mathbf{z})^2] &=& \frac{1}{n-1}  \sum_{i=1}^{n-1}12 + { 4\cos\left(2\frac{i}{n}\pi\right)} \overbrace{ + 16 \cos\left(\frac{i}{n}\pi\right)}^{\substack{\text{these terms cancel}\\\text{ due to  symmetry}}} \\
&=& \frac{1}{n-1} ( \sum_{i=1}^{n-1}12  + 4 \underbrace{\sum_{i=1}^{n-1}\cos\left(2\frac{i}{n}\pi\right))}_{=-1}  \\& =& \frac{12(n-1) -4}{n-1} \\&\approx& 12
\end{array}$$
where we used this to derive that the sum of cosines is equal to -1.

*We can not express the probability density function in a closed-form but we can express the cumulants of the distribution $\kappa_k(s(\mathbf{z}))$ in terms of the cumulants of a single chi-squared variable $\kappa_k(\chi_{(1)}^2)$. For increasing $n$ the 1st order cumulant will go to infinity (so to make the limiting distribution a normal distribution you should not only change the factor $1/(n-1)$ but also subtract the mean), the 2nd cumulant will approach to $12$ and the other, higher order cumulants, will approach zero, which means that you approach a normal distribution.
$$\kappa_k(s(\mathbf{z})) = \kappa_k(\chi_{(1)}^2)\frac{1}{\sqrt{n-1}^k} \sum_{i=1}^{n-1} \lambda_i^k \approx \kappa_k(\chi_{(1)}^2)\frac{n-1}{\sqrt{n-1}^k} \int_0^1  ( 2 + 2 \cos(x\pi))^k dx$$
Maybe there is a more direct way to use some version of the CLT for a sum of independent variables that only differ by a scaling constant instead of manually computing the cumulants. But I couldn't find one.
Why is it equivalent to a sum of chi-squared variables?
(see also here Show that the distribution of $x'Ax$ is linear combination of chi-squared)
Example for $n=3$
In the case of $n=3$ then the $z_1$ and $z_2$ are distributed like a multivariate normal distribution with a negative correlation. Geometrically it looks like an elongated shape.

n <- 10^4
set.seed(1)
x <- matrix(rnorm(3*n),3)
z1 <- x[2,]-x[1,]
z2 <- x[3,]-x[2,]
plot(z1,z2, xlab = expression(z[1]), ylab = expression(z[2]))

We can express the square in terms of alternative variables $Y_1 = \sqrt{0.5}(Z_1-Z_2) \sim N(0,3)$ and $Y_2 = \sqrt{0.5}(Z_1+Z_2) \sim N(0,1)$
$$Z_1^2 + Z_2^2 = 0.5(Z_1-Z_2)^2 + 0.5(Z_1+Z_2)^2  = Y_1^2 + Y_2^2$$
Note that the $Y_i$ are independent. So the distribution is similar to the distribution of a sum of independent squared normal distributed variables, but with different variance.
Generalized for all $n$
More generally $z_i$ is a multivariate normal distribution (any linear combination of the $z_i$ is a linear combination of the $x_i$ which is a normal distributed variable).
The variance of each $z_i$, being the sum of two standard normal variables, is $2$. The covariance of two neighbouring variables  is $-1$ (which you can find with covariance of sums). So the covariance matrix is like:
$$\Sigma = \begin{bmatrix} {2} & -1 & 0 & \dots & 0 &0\\
-1 & 2 & -1 & 0& \dots & 0  \\
0 & -1 & {2}& \dots & 0& 0\\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
0 & 0& 0 & \dots & {2} & -1 \\
0 & 0& 0 & \dots & -1 & {2}
\end{bmatrix}$$
In this general case we can do the same as for $n=3$ and restate the 
dependent $z_i^2$ as a sum of independent squared normal variables $y_i^2$. We use the same geometrical interpretation and rotate the distribution (keeping the radial distance invariant) and the distribution of $z_i$ is equivalent to rotated $y_i$ which have a variance that relate to the eigenvalues of the covariance matrix $\Sigma$. These eigenvalues will be in between 0 and 4 (see more about that below). 
These eigenvalues follow a cosine function
$$\lambda_i = 2 + 2 \cos(\frac{i}{n}\pi)$$
for $1\leq i\leq n+1$. Which can be derived from the general description of eigenvalues of triadagonal Toeplitz matrices (as mentioned by Aksakal in the comments, you can see previous edits of this post for an alternative derivation of that relation with cosines)
