This is basically the same idea as Nameless' answer, but I think it's cleaner if you just keep and update $SX$, the sum of the terms, and $SX^2$, the sum of their squares, update those via
$$\begin{align}
SX_{i+1} &= SX_i + x_{i+1}-x_{i-999}\\
SX^2_{i+1} &= SX^2_i + x_{i+1}^2-x_{i-999}^2\\
\end{align}
$$
and then compute
$$\begin{align}
mean_{i+1} &= SX_{i+1}/1000\\
var_{i+1} &= \frac{1}{1000}SX^2_{i+1} - (\frac{1}{1000}SX_{i+1})^2\\
&= \frac{1}{1000}SX^2_{i+1} - mean_{i+1}^2\\
\end{align}
$$
(This is the way that old calculators would let you accumulate an arbitrary number of data values and then let you compute mean and variance at the end.)
There is a possible issue in that you have to be able to store a number of magnitude roughly $1000*(\text{typical }x)^2$. Hopefully your values aren't too large.
If you're really worried about numerical instability, @awkward's link to Wellford's online algorithm looks to be the state-of-the-art; I'd check that out if your output values start to go wrong after you've been running for a while.