Fourth moments of Gaussian processes 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 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?
 A: 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. For many more details see  Janson's book Gaussian Hilbert Spaces
