Random walks in multinomial case

Question

Model:

A vector $X=(X_1, X_2, X_3)$ that follows a trinomial distribution with parameters $p=1/3$ and $n$.

(I have a coin with three sides $S1$, $S2$, $S3$). I flip the coin $n$ times. The coin has a probability $p=1/3$ to be flipped to the side $S1$, and similarly with $S2$ and $S3$. $X_1$ counts the number of times the coin is flipped to $S1$ ($X_2$ and $X_3$ are defined similarly).

Questions:

1. Let $0 \leq \alpha \leq n$. I Want to find $T(\alpha)= P (((X_1 - X_2) \geq \alpha) \cap ((X_1 - X_3) \geq \alpha))$, or a lower bound on $T(\alpha)$
2. For what values of $\alpha$, The lower bound on $T(\alpha)$ does not depend on $n$ (is a constant) ?

If $X$ was following a binomial distribution, the problem would have been easy to solve using random walks, but I do not know how to solve it in the multidimensional case.

Any idea? Thank you.

Two-dimensional case

$X= (X_1, X_2)$ follows a binomial distribution with parameters $p=1/2$ and $n$.

(A coin is flipped $n$ times. the coin is up with probability $p=1/2$ and down with probability $p$. $X_1$ counts the number of ups and $X_2$ counts the number of downs)

Let $Y= X_1 - X_2$ (We have $n$ coin flips, when the coin is flipped up, we add (+1) to $Y$, when the the coin is flipped down, we add (-1) to $Y$). We can see this as a random walk, when the coin is up we go to the right and when it is down we go to the left.

$\mathrm{Var}(Y)=n = \sqrt n ^2 $. $E(Y)=0$

By the Central Limit Theorem, $P(Y \geq \alpha) \approx 1 - \Phi(\alpha / \sqrt n)$ (normal distribution).

Answer 1

The probabilities for this problem can be calculated explicitly for quite large $n$.
To get a very good approximation for even moderately large $n$, we can use the multivariate central limit theorem.

Define $U_1 = X_1 - X_2$ and $U_2 = X_1 - X_3$. Note that by symmetry, $$ \newcommand{\e}{\mathbb{E}}\renewcommand{\Pr}{\mathbb{Pr}}\newcommand{\Cov}{\mathrm{Cov}}\e U_1 = \e U_2 = 0 \> . $$ We also have, by the bilinearity of the covariance operator, that $$ \Cov(U_1, U_2) = \Cov(X_1,X_1) - \Cov(X_1,X_3) - \Cov(X_1,X_2) + \Cov(X_2,X_3) \>. $$ Now, $\Cov(X_1,X_1) = n p (1-p)$ where $p = 1/3$ here. An only slightly more difficult calculation yields $$ \Cov(X_1,X_2) = - n p^2 \>, $$ and, of course, by symmetry again, $\Cov(X_1,X_2) = \Cov(X_1,X_3) = \Cov(X_2,X_3)$.

Hence, $\Cov(U_1, U_2) = n p (1-p) + n p^2 = np = n/3$ and by a similar calculation, $\Cov(U_1, U_1) = \Cov(U_2,U_2) = 2 n p = 2 n / 3$.

Observe that $U_1$ and $U_2$ are each the sums of independent and identically distributed random variables. For example, if $\xi_i \in \{1,2,3\}$ is the outcome of the $i$th draw, then $U_1 = \sum_{i=1}^n 1_{(\xi_i = 1)} - 1_{(\xi_i = 2)}$ where $1_{(\cdot)}$ is the indicator function.

Hence, by the multivariate central limit theorem, we conclude that $$ \sqrt{\frac{3}{n}} (U_1,U_2) \xrightarrow{d} \mathcal{N}(0, \Sigma) $$ where $$ \Sigma = \left(\begin{array}{cc}2 & 1 \\ 1 & 2\end{array}\right). $$

Now, since $$ T(\alpha) = \Pr( \{X_1 - X_2 \geq \alpha \} \cap \{X_1 - X_3 \geq \alpha \} ) = \Pr( U_1 \geq \alpha, U_2 \geq \alpha) $$ then, we can approximate $T(\alpha)$ as follows $$ T(\alpha) \approx \int_{\alpha\sqrt{3/n}}^\infty \int_{\alpha\sqrt{3/n}}^\infty \frac{1}{2 \sqrt{3} \pi} e^{-\frac{1}{3}(u_1^2 - u_1 u_2 + u_2^2 )} \mathrm{d}u_1 \mathrm{d}u_2 \> . $$

Below is some very brief $R$ code that compares a simulation of the true process against a simulation using the normal approximation assuming $n = 100$ underlying trinomial trials.
First, the picture.

Here is the code.

set.seed(.Random.seed[1])
n <- 100
N <- 10000

X   <- matrix( sample(1:3, n*N, replace=T), nc=n )
xt  <- apply(X,1,table)
dxt <- cbind( xt[1,]-xt[2,], xt[1,]-xt[3,] )
xtt <- table(apply(dxt,1,min))

L   <- matrix( c(sqrt(2), 1/sqrt(2), 0, sqrt(3/2)), 2 )
Y   <- L %*% matrix( rnorm( 2*N ), nr=2 ) * sqrt(n) / sqrt(3)
ytt <- table( floor(apply(Y,2,min)) )

plot( names(xtt), xtt/N, type="h", xlab="a", ylab="Density of T(a)" )
lines( names(ytt), ytt/N, col="red", type="h" )
legend( "topright", legend=c("actual", "normal approx."), lty="solid",
        col=c("black", "red"), bg="white", inset=0.02 )

There are also R packages to calculate bivariate normal densities and probabilities. Both mnormt and fMultivar are examples, but I don't know enough to recommend any one over another.