Calculate the expectation of the product of two random variables There is a sequence of independent and identically distributed continuous random variables $X_1,X_2,\ldots$ with common density function $f(x)$. We say that a record occurs at time $n$ if $X_n > \max(X_1,X_2,\ldots,X_{n-1})$. Consider the random variable $Y_i$ defined as: $Y_i = 1$ if a record occurs at time $i$, $Y_i = 0$ otherwise. Compute $E[Y_i \cdot Y_j]$ where $i < j$.
I compute like this:
$$E[Y_i \cdot Y_j] = P(Y_i = Y_j = 1) = P(Y_i=1 \mid Y_j=1)P(Y_j=1).$$
Since $P(Y_j=1)=1/j$ and $$P(Y_i=1 \mid Y_j=1) = P(Y_i=1) = 1/i,$$ then $E[Y_i \cdot Y_j]=1/(ij)$.
It seems not quite right and I am not very sure about it. The reason why I'm not so sure about it is that, it seems by this way, the events "the record occurs at the time $i$" and "the record occurs at the time $j$" would be independent. But given "the record occurs at the time $j$", would it be harder for "the record occurs at the time $i$" to happen? Because $X_j > X_i$, and $X_i$ cannot be very big to be the record.
Thank you for your help!
 A: You should be able to express the probability as (see below)
$P(Y_i = Y_j = 1) = \int_{(0,1)^2} p(u_1,u_2) \mathrm{d}(u_1,u_2)$,
where for a iid sequence $(U_k)$, each $\text{Unif}(0,1)$-distributed,
$p(u_1,u_2)= P(u_2 \geq u_1, u_2 \geq U_{j-1}, \ldots , u_2 \geq U_{i+1}, u_1 \geq U_{i-1}, \ldots , u_1 \geq U_1)$.
This probability is due to independence,
$p(u_1,u_2) = \mathbb{1}_{u_2\geq u_1} \prod_{k=i+1}^{j-1} P(U_k\leq u_2) \prod_{k=1}^{i-1} P(U_k \leq u_1)$,
where the product as usually is $1$ if the lower index is higher than the upper index. Then the integral becomes
$P(Y_i = Y_j = 1)
=
\int_{(0,1)} u_2^{j-i-1}  \int_0^{u_2} u_1^{i-1}    \, \mathrm{d}u_1 \mathrm{d}u_2, \\
=
\int_{(0,1)} \frac{1}{i} u_2^{j-1} \, \mathrm{d}u_2
=
\frac{1}{ij},
$
in accordance with your intuition. I am actually quite surprised by this result - maybe there is an easier way to show it. But your argument of independence is in my opinion not clear enough for a formal proof. Also, the above does not imply independence.
A: I tried to run some simulations to figure out, what was going on. The interested reader could try to run the following R-code, where I have considered some different scenarios.
rm(list = ls())

set.seed(2)

### Uniform case ###

n <- 20
m <- 10000

U <- matrix(runif(n*m), nrow = m, ncol = n)

Y <- matrix(0, nrow = m, ncol = n)

for(k in 1:m){
  Y[k, ] <- c(1, (U[k, 2:n]>cummax(U[k, 1:(n-1)])))
}

#Some comparisons:
mean(Y[ ,5]*Y[ ,7]) ; 1/(5*7)
mean(Y[ ,2]*Y[ ,13]); 1/(2*13)
mean(Y[ ,3]*Y[ ,4]) ; 1/(3*4)
mean(Y[ ,16]*Y[ ,18]) ; 1/(16*18)

### Non-uniform case ###

n <- 20
m <- 10000

E <- matrix(rexp(n*m, 2), nrow = m, ncol = n)

Y <- matrix(0, nrow = m, ncol = n)

for(k in 1:m){
  Y[k, ] <- c(1, (E[k, 2:n]>cummax(E[k, 1:(n-1)])))
}

#Some comparisons:
mean(Y[ ,5]*Y[ ,7]) ; 1/(5*7)
mean(Y[ ,2]*Y[ ,13]); 1/(2*13)
mean(Y[ ,3]*Y[ ,4]) ; 1/(3*4)
mean(Y[ ,16]*Y[ ,18]) ; 1/(16*18)

### Discrete case ###

n <- 20
m <- 100000

D <- matrix(rbinom(n*m, 50, 1/2), nrow = m, ncol = n)

Y <- matrix(0, nrow = m, ncol = n)

for(k in 1:m){
  Y[k, ] <- c(1, (D[k, 2:n]>cummax(D[k, 1:(n-1)])))
}

#Some comparisons:
mean(Y[ ,5]*Y[ ,7]) ; 1/(5*7)
mean(Y[ ,2]*Y[ ,13]); 1/(2*13)
mean(Y[ ,3]*Y[ ,4]) ; 1/(3*4)
mean(Y[ ,16]*Y[ ,18]) ; 1/(16*18)

