Independent Bernoulli trials vs markov chain 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.

*

*Consider the trials, $X_t \stackrel{iid}{\sim} Bern(\theta)$ and estimate $\theta$.

*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.
Partial Answer
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?
 A: 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.
