There is a section entitled 'Bivariate Half-normal distribution in:
Continuous Multivariate Distributions: Models and applications
By Samuel Kotz, Norman Lloyd Johnson, N. Balakrishnan.
I would be curious to see how this can be generalized to a random vector of any dimensions.
In fact, the bivariate case appears to be thoroughly treated in this paper:
http://www.stat-athens.aueb.gr/~jpan/papers/Panaretos-ApplStatScience2001(119-136)ft.pdf
Suppose $A=((a_{ij}))$ and $y=(y_1,y_2,\ldots,y_n)'\sim N(0,I_n)$, so that $y_i$s are i.i.d standard normal.
You already have $\operatorname E[y'Ay]=\operatorname{tr}(A)$ by the result you quoted (discussed here).
For the variance, you can simply use $$\operatorname{Var}(y'Ay)=\operatorname E[(y'Ay)^2]-(\operatorname E[y'Ay])^2\tag{1}$$
To compute the first expectation, we write the quadratic form as $y'Ay=\sum_{i,j}a_{ij}y_iy_j$ to get $$(y'Ay)^2=\sum_{i,j,k,l}a_{ij}a_{kl}y_iy_jy_ky_l \tag{2}$$
Now observe that $$\operatorname E[y_iy_jy_ky_l]=\begin{cases}3&,\text{ if }i=j=k=l \\ 1&,\text{ if }i=j,k=l;i=k,j=l;i=l;j=k \\ 0&,\text{ otherwise }\end{cases}$$
Therefore taking expectation on both sides of $(2)$,
$$\operatorname E[(y'Ay)^2]=3\sum_i a_{ii}^2+\sum_i\left(\sum_{k\ne i}a_{ii}a_{kk}+\sum_{j\ne i}a_{ij}^2+\sum_{j\ne i}a_{ij}a_{ji}\right)\tag{3}$$
Keeping in mind that $A$ is symmetric, you have $\operatorname{tr}(A^2)=\sum_{i,j}a_{ij}^2$.
It is now straightforward to see that $(3)$ reduces to $$\operatorname E[(y'Ay)^2]=(\operatorname{tr}(A))^2+2\operatorname{tr}(A^2)$$
From $(1)$ you get the desired result $$\boxed{\operatorname{Var}(y'Ay)=2\operatorname{tr}(A^2)}$$
You can now try to generalize this to the case $y\sim N(0,\Sigma)$ and hence to $y\sim N(\mu,\Sigma)$.
Reference:
- Linear Regression Analysis by Seber and Lee.
Best Answer
The variance is the matrix product: $$ 1'\Sigma1 $$