# Original Question

Suppose we have a sequence of Bernoulli trials $$X_1, X_2, \cdots X_T$$ which are ordered in time and may or may not be independent. I am interested in understanding the probability of success.. The way I'm thinking about it, I have two primary options.

1. Consider the trials, $$X_t \stackrel{iid}{\sim} Bern(\theta)$$ and estimate $$\theta$$.
2. Treat the observations as a markov chain (does success/failure at time $$t$$ have an effect on probability of success at time $$t+p$$?

How do i determine which model is best? Or in otherwords, how do I estimate the order of the markov chain.

To test the validity of the first model, I was considering using the approach of Andrew Gelman in his Bayesian Data Analysis book, where he constructs an "empirical p-value" via simulation using the test statistic $$T = \text{longest streak of successes}$$. Beyond this, I am unsure how to appropriately compare the different models.

It appears that the PACF (partial autocorrelation function) can be used to estimate the order of these sequences.

I generated 3 large sequences $$(T=1000)$$ of Bernoulli trials, under independence, first order Markov dependence and second order Markov dependence. We plot the ACF for each case.

Moreover, for the first order Markov Chain, the ACF seems to satisfy the property of an AR(1) model that $$\rho_k = \rho_1^k$$.

I can't seem to find any formal justification for why this works however. As I understand it, if the Bernoulli trials satisfy a $$p$$ order Markov property, this does not imply the chain is an $$AR(p)$$ process, ecause we cannot write $$X_t = \phi X_{t-1} + \epsilon_t$$ for $$\epsilon_t$$ iid. (If somebody can confirm or correct me on this, that would be great).

So if the $$p$$ order Markov process is not an $$AR(p)$$ process, why does this work?

• First you say that $X_1,...,X_T$ may be dependent and then you assume they are i.i.d. -- could you clarify on that? – Tim Oct 27 '17 at 20:33
• Sorry for the confusion. I'm not assuming they're iid. I'm saying that since I can't write the process as $X_t = \phi X_{t-1} + \epsilon_t$ for $\epsilon_t$ iid, then it's not an AR(1) process (I think). – knrumsey Oct 27 '17 at 20:38
• Unless you're talking about the very beginning. I say that they may or may not be independent. So I definitely want to explore the iid case, but I can't rule out a first or second order markov model. – knrumsey Oct 27 '17 at 20:47
• In the first case you have a simple beta/binomial model. That is if you assume $X_i$ are independent conditioned on $\theta$. – Łukasz Grad Oct 27 '17 at 21:36
• Yes, I know how to handle each case. The question is about determining an appropriate model, not parameter estimation. – knrumsey Oct 27 '17 at 21:42

As a general rule, if you want to test a hypothesis $H_0$ against $H_A$, you must specify a model that encapsulates both possibilities and then use that model to try to infer which of the hypotheses is correct. So if you think your observed time-series data might have some form of auto-correlation, and you want to test this, it is no good modelling them as IID. You need to use a model that allows for auto-correlation, but also allows independence as a special case, so that these competing hypotheses can be tested from within the model.

You are correct that the standard AR(1) auto-correlation structure for a continuous variable is no good in this case. You will instead need to formulate some kind of model that is suitable for a sequence of Bernoulli trials. For a sequence of Bernoulli trials, a stationary Markov chain for this process would have two parameters $0< \theta_0 <1$ and $0<\theta_1<1$ and use the recursive equations:

$$\begin{matrix} X_1 \sim \text{Bern} \left( \frac{\theta_0}{1+\theta_0-\theta_1} \right) & & & X_{t+1} | X_t \sim \text{Bern}(\theta_{X_t}). \end{matrix}$$

This gives a general stationary Markov chain, where the probability in the Bernoulli trial depends on the previous outcome. It can be shown that $\mathbb{Corr}(X_{t+1}, X_t) = \theta_1 - \theta_0$ within this model, so if you want to test for independence, you would be testing the hypotheses:

$$\begin{matrix} H_0: \theta_1 = \theta_0 & & H_A: \theta_1 \neq \theta_0. \end{matrix}$$

It should be possible to fit this model to your data, estimate the parameters of the model, and test the hypothesis of independence. This could be done using classical or Bayesian methods.

• Thank you! Can I just ask what motivates the success probability for $X_1$? Why not an average (or weighted average) of $\theta_1$ and $\theta_2$. Under $H_0$, this would give the same model, but seems more interpretable to me, unless I'm missing something. – knrumsey Jan 29 '18 at 19:18
• Great question. The expression you're looking at is the marginal probability of success that comes from the recursive equation under the condition of stationarity. To see that, let $\phi \equiv \mathbb{P}(X_t = 1)$ for all $t$. Apply the law of total probability using the recursive equation to get $\phi = \phi \theta_1 + (1-\phi) \theta_0$. Now solve for $\phi$. – Ben Jan 29 '18 at 21:13
• If you give a different marginal success probability for $X_1$ then the recursive equation will operate to change the marginal success probability over time, and so you won't have a stationary model. – Ben Jan 29 '18 at 21:15
• thank you again. This model worked like a charm with some beta priors for $\theta_0$ and $\theta_1$. How would you go about showing that the correlation is $\theta_1 - \theta_0$. Secondly, do you have any thoughts on how to extend this model for Markov chains of higher order? – knrumsey Feb 23 '18 at 21:26
• I'm guessing the equivalent approach (for second order Markov) would be to define $\theta = (\theta_{00}, \theta_{01}, \theta_{10}, \theta_{11})$ and define $X_t \sim Bern(\theta_{X_{t-1}X_{t-2}})$ for $t \geq 3$. But how do we define the first two success probabilities in order to ensure stationarity? – knrumsey Feb 23 '18 at 21:59