Variance of the Product of Correlated Random Variables?

Question

I would like to multiply two correlated random variables, but I'm getting a negative variance. Please point out where I'm wrong.

Variable1 and Variable2 are projected onto data2 from a model trained on data1. The equation I am trying to use to calculate the variance of their products is from Variance of product of dependent variables. But when I translate it into R I get:

> data2$ProductVariance <- cov(data1$Variable2^2, data1$Variable1^2, use = "na.or.complete") + (data2$Variable2.se^2 + data2$Variable2^2)*(data2$Variable1.se^2 + data2$Variable1^2) - (cov(data1$Variable2, data1$Variable1, use = "na.or.complete") + data2$Variable2*data2$Variable1)^2
> data2[c(1, 50, 100, 150, 200), c("Variable1", "Variable1.se", "Variable2", "Variable2.se", "ProductVariance")]
    Variable1 Variable1.se  Variable2 Variable2.se ProductVariance
2   102.32145    0.4053362 0.68919721  0.006114099       -3.649645
60  103.55957    0.5298381 0.66190120  0.006577132       -3.528850
120  99.81072    0.4176735 0.07347951  0.005274879       -3.735716
183 100.59532    0.4008085 0.62787019  0.005122185       -3.777981
246 102.27556    0.3762328 0.73455316  0.006578060       -3.597734

Their covariance and squared-covariance are:

> cov(data1$Variable2, data1$Variable1, use = "na.or.complete")
[1] 0.0008326011
> cov(data1$Variable2^2, data1$Variable1^2, use = "na.or.complete")
[1] -4.00164

And, broken out by term, it looks like the negatives are coming from the 1st and 3rd term.

> data2$FirstTerm <- cov(data1$Variable2^2, data1$Variable1^2, use = "na.or.complete")
> data2$SecondTerm <- (data2$Variable2.se^2 + data2$Variable2^2)*(data2$Variable1.se^2 + data2$Variable1^2)
> data2$ThirdTerm <-  -(cov(data1$Variable2, data1$Variable1, use = "na.or.complete") + data2$Variable2*data2$Variable1)^2
> data2[c(1, 50, 100, 150, 200), c("Variable1", "Variable1.se", "Variable2", "Variable2.se", "ProductVariance", "FirstTerm", "SecondTerm", "ThirdTerm")]
    Variable1 Variable1.se  Variable2 Variable2.se ProductVariance FirstTerm SecondTerm  ThirdTerm
2   102.32145    0.4053362 0.68919721  0.006114099       -3.649645  -4.00164 4973.49137 -4973.1394
60  103.55957    0.5298381 0.66190120  0.006577132       -3.528850  -4.00164 4699.16893 -4698.6961
120  99.81072    0.4176735 0.07347951  0.005274879       -3.735716  -4.00164   54.06632   -53.8004
183 100.59532    0.4008085 0.62787019  0.005122185       -3.777981  -4.00164 3989.61583 -3989.3922
246 102.27556    0.3762328 0.73455316  0.006578060       -3.597734  -4.00164 5644.57031 -5644.1664

Have I misunderstood the formula?

EDIT

Should

> cov(data1$Variable2^2, data1$Variable1^2, use = "na.or.complete")
[1] -4.00164

be positive?

Answer 1

Expanding on my earlier comment, here's the result of just switching out your uses of $x_i$ and $y_i$ for $E(X)$ and $E(Y)$:

# Set up x and y
x <- c(102.32145, 103.55957, 99.81072, 100.59532, 102.27556)
y <- c(0.68919721, 0.66190120, 0.07347951, 0.62787019, 0.73455316)
# Put them into a dataframe and take their squares
df <- data.frame(x, y)
df$x.sqr <- (df$x)^2
df$y.sqr <- (df$y)^2
# Find the variances of x and y
x.var <- var(df$x)
y.var <- var(df$y)
# Find the expected values of x and y
x.exp <- mean(df$x)
y.exp <- mean(df$y)

# 1st term, find cov(x.sqr, y.sqr)
first.term <- cov(df$x.sqr, df$y.sqr)
# first.term = 48.1506

# 2nd term, find [V(x) + E(x)^2]*[V(y) + E(y)^2]
a <- x.var + (x.exp)^2
b <- y.var + (y.exp)^2
second.term <- a*b
# second.term = 3987.9941

# 3rd term, find [cov(x,y) + E(x)E(y)]^2
c <- cov(df$x, df$y)
d <- x.exp * y.exp
third.term <- (c + d)^2
# third.term = 3248.7993

# All together now!
var.product <- first.term + second.term - third.term
# var.product = 787.3454

And the variance is positive! Obviously in this case the variance is very high - I've calculated it using just the five rows of data you display above (i.e. $n = 5$). I would imagine when you compute it for your whole sample size with a much larger value of $n$ the variance will be much more reasonable.