I am currently trying to deduce the theoretical distribution of the stochastic integral $$ \int_0^1 2W_t \hspace{1mm} dW_t $$
To attempt this, I have first evaluated the above to give $$ (W_1)^2 - 1 \hspace{10mm} (*) $$
Next, I know that $W_1 \sim N(0,1)$, and thus $(*)$ would also be normally distributed with mean $$ \mathbb{E} [W_1^2 - 1] = \mathbb{E} [W_1^2] - 1 = 0 $$ and variance $$ \text{Var}(W_1^2 - 1) = \text{Var}(W_1^2) = \mathbb{E} [W_1^4] - \mathbb{E} [W_1^2]^2 = 3-1 = 2 $$
However, I also have a script in R which simulates 1000 Wiener processes and evaluates the integral above based on each of these. The result of this is an array containing 1000 values, which have mean $-0.03$ and variance $0.44$. I am thus wondering whether I have calculated the theoretical variance incorrectly?
Your derivation
$$ \int_0^t 2W_t dW_t = (W_1)^2 -1$$
is correct and follows from a simple application of Ito's formula.
It follows that the expectation and variance are given by the expressions you give
$$\mathbf E \left[ \int_0^1 2W_t dW_t \right] = 0, \qquad \text{Var}\left( \int_0^1 2W_t dW_t \right) = 2.$$
The question then is how you are simulating the stochastic integral $\int_0^1 2W_t dW_t$. Two immediate thoughts come to mind:
At how many breakpoints are you computing the integral?
Assuming that you are approximating the integral by computing
$$I_n = 2\sum_{k=0}^{n-1}\textstyle W\left(\frac{k}{n}\right) \left( W\left(\frac{k+1}{n}\right) - W\left(\frac{k}{n}\right) \right),$$ then this will have a calculation error associated with it due to the approximation; in fact you can calculate
\begin{align*} \text{Var} \left( I_n\right) & = 4 \sum_{k=0}^{n-1} \text{Var} \bigg( \textstyle W\left(\frac{k}{n}\right) \left( W\left(\frac{k+1}{n}\right) - W\left(\frac{k}{n}\right) \right) \bigg) \\ & = 4 \sum_{k=0}^n \text{Var}\textstyle \big( N\big(0,\frac{k}{n}\big)\text{Var}\textstyle \big( N\big(0,\frac{1}{n}\big) \big) \\ & = 4 \sum_{k=0}^n \frac{k}{n^2} \\ & = 2 - \frac2n. \end{align*}
Rearranging this, your calculation of the variance being 0.44 would suggest you are using $n=1,2$ which seems unlikely!
Are you sure you are not in fact computing $\int_0^1 2W_t dt$.
If this were the case, then similar stochastic calculus implies
$$\mathbf E \left[ \int_0^1 2W_t dt \right] = 0, \qquad \text{Var}\left( \int_0^1 2W_t dt \right) = \frac43.$$
When implementing this computationally, further approximation errors will be induced as described for the previous method.
Perhaps could you explain in more detail how you are doing the computation?