I am working on a topic outside my main research area, so I am afraid I am reproving obvious results, so I would like to ask for a reference. Google didn't help, mostly because I am looking for formulas.


Let $v$ be an $n$-variate Gaussian random variable, with $E[v]=0$ and $E[vv^\top]=I$. If I am not mistaken, it follows (with a bit of work) from [Isserlis' theorem](https://en.wikipedia.org/wiki/Isserlis%27_theorem) that for any $a,b\in\mathbb{R}^n$
$$
E[(a^\top v)^2(b^\top v)^2] = (a^\top a)(b^\top b) + 2 (a^\top b)^2.
$$
**Q1. Do you have a reference for this identity?**


Now let us consider a generalization: $v\_t$, for $t\in \mathbb{Z}$, is a $n$-variate random process, i.e., each $v\_t$ is a vector-valued random variable, and $E[v\_t]=0$, $E[v\_tv\_t^\top]=I$, $E[v\_tv\_s^\top]=0$ for $s\neq t$. Let $L$ be a shift operator on the sequence, i.e., $Lv\_t = v\_{t-1}$, and $a(L)$ be a polynomial in it with coefficients $a\_i\in\mathbb{R}^n$, i.e., $a(L)v\_t = a\_0 v\_t + a\_1v\_{t-1} + \dots + a\_d v\_{t-d}$. Then, I imagine that there should be a result of the kind
$$E[(a(L)^\top v)^2(b(L)^\top v)^2] = \text{the coefficient of $L^0$ in } (a(L)^\top a(L^{-1}))(b(L)^\top b(L^{-1})) + 2 (a(L)^\top b(L^{-1}))^2.$$


**Q1. Is a formula of this kind true / known? Where can I find references on this topic and this kind of algebra on processes?**



Fourth moments of Gaussian processes
The first equality you mention is a special case of *Wick's formula* or *diagram formula*. Suppose that you have a Gaussian random vector $X=(X\_1,\dotsc, X\_n)$ that is *centered*, i.e., $\newcommand{\bE}{\mathbb{E}}$ $\newcommand{\bR}{\mathbb{R}}$


$$ \bE[X\_i]=0,\;\;\forall i=1,\dotsc, n. $$


A special case of Wick's formula computes $\bE[X\_1\dotsc X\_n]$. Note that this expectation is $0$ if $n$ is odd so we only need to consider the case $n=2k$.


A *Feynman graph* associated to $\{X\_1,\dotsc, X\_{2k}\}$ is then a graph with vertex set $\{X\_1,\dotsc, X\_{2k}\}$ and exactly $k$-edges so no two edges have a vertex in common. If $\newcommand{\be}{\boldsymbol{e}}$ $\be$ is an edge of a Feynman graph with endpoints $X\_i, X\_j$ then we define its *weight* to be the covariance $\bE[X\_iX\_j]$. The weight of a Feynman graph $\Gamma$ is then


$$w(\Gamma):=\prod\_{\be \in E(\Gamma)}w(\be), $$


where $E(\Gamma)$ denotes the set of edges of $\Gamma$. Wick's formula states that


$$ \bE[X\_1X\_2\dotsc X\_{2k}]=\sum\_{\Gamma} w(\Gamma) $$


where the summation is carried over the set of Feynman graphs with vertex set $\{X\_1,\dotsc, X\_{2k}\}$.


Here is how it applies to your case. Suppose that $Y=(Y\_1,Y\_2)$ is a centered Gaussian vector. Then the random vector 


$$ X=(X\_1,X\_2,X\_3,X\_4)=(Y\_1,Y\_1,Y\_2,Y\_2) $$


is also a centered random Gaussian vector and


$$ \bE[X\_1X\_2X\_3X\_4]=\bE[Y\_1^2Y\_2^2]. $$


In this case there are $3$ Feyman graphs:


$$ \Gamma\_1,\;\;E(\Gamma\_1)= \big\{ (1,2), (3,4)\big\},\;\; w(\Gamma\_1)= \bE[Y\_1^2]\bE[Y\_2^2], $$


$$\Gamma\_2,\;\;E(\Gamma\_2)=\big\{ (1,3), (2,4)\big\} ,\;\; w(\Gamma\_2)=\bE[Y\_1Y\_2]^2, $$


$$\Gamma\_2,\;\;E(\Gamma\_2)=\big\{ (1,3), (2,3)\big\} ,\;\; w(\Gamma\_3)=w(\Gamma\_2)=\bE[Y\_1Y\_2]^2. $$


Wick's formula then implies


$$ \bE[Y\_1^2Y\_2^2]= \bE[Y\_1^2]\bE[Y\_2^2]+2\bE[Y\_1Y\_2]^2. $$


The result you mention in $Q\_1$ corresponds to $Y\_1= a^\top v$, $Y\_2=b^\top v$.


The fact that $(Y\_1,Y\_2)$ is a Gaussian vector follows from the following property of Gaussian vectors: if $V=(V\_1,\dotsc, V\_n)$ is a Gaussian vector and $A: \bR^n\to\bR^m$ is a linear map, then $AV$ is an $m$-dimensional Gaussian vector. 


If we choose $A: \bR^n\to \bR^2$, $V\mapsto (Y\_1,Y\_2)=(a^\top V, b^\top V)$ then you see that $(Y\_1,Y\_2)$ is Gaussian.


For the generalization you mention, there are plenty to say if the process is Gaussian.


For details see Sections 1 and 2 of [these notes](http://www3.nd.edu/~lnicolae/Ito_integral.pdf). For many more details see Janson's book *Gaussian Hilbert Spaces*


