I have two realizations of a poisson stochastic process, they are over the same space with rate $\lambda_{1}$ and $\lambda_{2}$. What is the probability that N elements in both sequences are the same, e.g., $X(t) \in [0,1000]$?
Note: The two processes are independent, i.e. X(s) is independent of X(t). Lets say the number of elements within a process is finite, e.g. K.
The question describes two independent sequences of $K$ iid Poisson variates. Let's use indexes $1,2$ to distinguish them, so that $X_1$ has parameter $\lambda_1$ and $X_2$ has parameter $\lambda_2$.
Because the question asks for probabilities, I presume we do not yet have any realizations: otherwise the answer to any question of "probability" is either $0$ or $1$, since the assertion is either false or true, respectively. Instead, it is desired to know beforehand what the chances are that the two processes will agree [exactly] $N$ times.
This can easily be worked out from properties of the Skellam distribution, which is a difference of two independent Poisson processes: we just compute the chance that it equals $0$ and proceed from there. However, it's often useful to solve problems of probability from first principles, so let's proceed as if we did not know of the Skellam distribution.
At each time $t$, $t = 1, 2, \ldots, K$, the chance that $X_1(t)=X_2(t)$ is the chance that they both equal $0$ plus the chance they both equal $1$ plus ... etc, whence--summing over this exhaustive partition of the event--it is immediate from the Poisson distribution formula that
$$\Pr(X_1(t) = X_2(t)) = \exp(-\lambda_1 - \lambda_2) \sum_{i=0}^\infty \frac{(\lambda_1 \lambda_2^i)^i}{i! i!} = e^{-\lambda_1-\lambda_2} I_0\left(2 \sqrt{\lambda_1\lambda_2}\right).$$
($I_0$ is a modified Bessel function of the first kind of order zero.)
This figure shows contours (on a logarithmic scale) of the chance that $X_1(t)=X_2(t)$. The axes are $\lambda_1$ and $\lambda_2$.
The distribution of the number of such equalities over $K$ independent trials has (by definition) a Binomial distribution with parameters $p = e^{-\lambda_1-\lambda_2} I_0\left(2 \sqrt{\lambda_1\lambda_2}\right)$ and $K$, whence the chance of exactly $N$ equalities equals $\binom{K}{N}p^N(1-p)^{K-N}$. The chance of $N$ or more equalities is similarly computed from the complementary CDF of the Binomial distribution (given by a beta integral) or, more simply, by summing the chances for $N$, $N+1$, ..., up to $K$.
If the example in the question is intended to assert that the variables are truncated to the range $[0, 1000]$, then (1) the summation should extend only from $i=0$ to $i=1000$ and (2) the normalizing factors $\exp(-\lambda_i)$ must be replaced by truncated series $\sum_{j=0}^{1000} \lambda_i^j / j!$. Generalizing the value $1000$ to, say, $n$, we have
$$\Pr(X_1(t) = X_2(t)) = e^{-\lambda_1 - \lambda_2}C_1 C_2\left(I_0\left(2 \sqrt{\lambda }\right)-\frac{\lambda ^{n+1} \, _1F_2(1;n+2,n+2;\lambda )}{((n+1)!)^2}\right)$$
where $C_i = \frac{\Gamma (n+1,\lambda_i )}{\Gamma (n+1,\lambda_i )-\Gamma (n+1)}$, $\lambda = \sqrt{\lambda_1 \lambda_2}$ and $_1F_2$ is a generalized Hypergeometric function. $\Gamma (n+1,\lambda_i)$ are incomplete Gamma functions (where the integration begins at $\lambda_i$ rather than at $0$). Expressing the truncated sums in this form provides access to alternative (possibly faster, more accurate) computational methods compared to actually computing the sums.
