Let $X$ and $Y$ be independent random variables each having the same geometric distribution. We need to find $P(X=Y)$.
I tried as follows :
Let us say $X$~$Geo(p)$ and $Y$~$Geo(\theta)$.
$P(X=Y) = \sum_{y} P(X=Y|Y=y)P(Y=y) =>\sum_{y} P(X=y)P(Y=y)$.
=> $\sum_{y} p(1-p)^{y-1} \theta(1- \theta)^{y-1}$ => $\sum_{y} (p \theta)[(1-p)(1- \theta)]^{y-1} $ => $ \dfrac{p \theta}{1-(1-p)(1- \theta)}$
Is this correct ? (The independence condition is never used ).
First, there are two styles of geometric random variable. Yours counts the trials until the first Success, the other counts Failures until the first Success.
Second, 'iid' means independent and identically distributed, so that you need to have $p = \theta.$ In case $p = \theta = 1/3,$ your answer computes to $1/(9-4) = 1/5 = 0.200.$
Finally, from the discussion in the Comments, I am not entirely certain you believe you are on the right track. So I will show you results of a simulation in R statistical software with a million realizations of $X$ and $Y$ with $p = 1/3.$
Notes: (1) Because R uses the 'other' version of the geometric random variable, I add $1$ to each realization in order to match the model you are using. (2) The comparison vector
x == yhas a million entries, each of then eitherTRUE(for equality) orFALSE. Themeanof this vector is its proportion ofTRUEs. (3) With a million iterations, the answer should be correct to three decimal places. So the answer is consistent with your derivation.