7 Further Topics in Probability

There are two important topics in probability that we haven’t yet managed to cover. One is stochastic process theory and the other is simulation. Let’s go over some of the key results.

7.1Stochastic Processes¶

A stochastic process on $\RR^K$ is a sequence of random vectors $\{\boldx_t\}_{t\geq 0}$ defined on a common probability space $(\Omega, \fF, \PP)$. The definition has relatively little content, but the fact that all these random vectors live on the same probability space should not be overlooked. It means that probabilities such as $\PP\{\boldx_t \geq 0 \text{ for all } t\}$ or $\PP\{\boldx_t \leq -10 \text{ for some } t\}$ are well defined. Putting the elements of the process on a common probability space pins down the joint distribution across $t$.

To do any kind of analysis with stochastic processes, we need to put some structure on the joint distribution, or, equivalently, on the random elements. In practice, this is almost always done recursively. Examples are given below. For the moment, however, let’s look at stochastic processes in an abstract setting and state some definitions.

One kind of stochastic process is an IID sequence. Such sequences are ideal to work with in a statistical setting. For example, the law of large numbers tells us that in an IID environment, finite expectations can always be approximated with sample averages. If we want to know how good the approximation is, we can look to the central limit theorem.

Once we move out of the IID world, however, these helpful properties can break down. Without restrictions on our stochastic process, we can say very little about estimation. The first task, then is to describe some properties that make our stochastic process resemble an IID process in some way.

7.1.1Stationarity and Ergodicity¶

As a first step, let’s define a type of process that preserves one of the properties of IID processes (the "ID" part). A stochastic process $\{ \boldx_t \}_{t \geq 0}$ is called stationary if it has the property that the distribution of any subset of these random vectors is unaffected by shifting them forward in time. That is,

for any $m \in \NN$ and any sequence of integers $t_1, \ldots, t_k$.

Clearly, any stationary sequence is identically distributed. In other words, for any stationary sequence, the sequence of marginals $\{\lL(\boldx_t)\}$ is a constant sequence.

Stationarity alone isn’t strong enough to give us properties like the LLN and CLT.

The problem with the process in example 7.1.4 is excessive dependence. New realizations provide insufficient (in fact no) new information. On the flip side, dependent processes can satisfy versions of the LLN if dependence is not too excessive (e.g., if $x_t$ and $x_{t+k}$ are "almost independent" for $k$ sufficiently large.) These are the so-called ergodic processes.

The formal definition of ergodicity varies across different literatures. Here we define a stationary stochastic process $\{\boldx_t\}$ on $\RR^K$ to be ergodic if the LLN holds; that is, if

for any $\bB$-measurable $h \colon \RR^N \to \RR$ such that $\EE |h(\boldx_t)| < \infty$. (Since $\{\boldx_t\}$ is stationary, this expectation does not depend on $t$.) One of the main tasks of this chapter will be deriving sufficient conditions for ergodicity in the kinds of processes economists often work with.

Incidentally, a casual definition of ergodicity one often hears is that "cross sectional and time series averages coincide." What does that mean in the present context?

To answer the question, imagine that $\boldx_t = x_t$ is a binary variable with $x_t = 1$ if a given individual is employed at time $t$ and zero otherwise. In this case, $\frac{1}{T} \sum_{t=1}^T x_t$ is the fraction of time spent employed over the period $1, \ldots, T$. With $h(x)=x$, the right-hand side of (7.2) is $\EE x_t = \PP\{x_t = 1\}$, or the stationary probability of employment. If we observe time $t$ employment outcomes $x_t^1, \ldots, x_t^N$ of individuals $1, \ldots, N$ and, moreover, the individuals all follow the same model and are sufficiently independent, the cross sectional average $\frac{1}{N} \sum_{n=1}^N x_t^n$ will also be close to $\PP\{x_t = 1\}$. Hence, with ergodicity and large samples, cross-sectional and time series averages (almost) coincide.

7.1.2Stochastic Recursive Sequences¶

In this section let’s turn to a more concrete view of stochastic processes. Afterward we’ll circle back and consider issues like ergodicity again.

Figure 7.1 shows the evolution of US GDP per capita since the early postwar era. The top panel is per capita GDP itself, while the second panel is the growth rate. How could we go about modeling such time series in a parsimonious way?

No reasonable model can be deterministic, since shocks play an important role in year-on-year change in GDP, and an IID process is too simple. We can, however, get close to both of these time series with one of the most common nontrivial stochastic processes: the scalar Gaussian AR(1) model

Here $a, b$, and $c$ are parameters, while $x_t$ is called the state variable. Equation (7.3) is an example of a stochastic difference equation. The process $\{x_t\}$ it defines is called a stochastic process or stochastic recursive sequence. A realization of the process is called a time series.

The dynamics of $\{x_t\}$ depends on the parameters. We can get some sense of this from simulation. Six individual time series $\{x_t\}$ are shown in figure 7.2, each generated using a different value of $a$. The values of $b$ and $c$ are both held fixed at 1, while $x_0 = 0$. R code for generating the figure is given in listing 3.

As suggested by the figure, the simulated time paths are sensitive to the value of $a$. If $a$ is outside the interval $(-1,1)$, the series tend to diverge. If $|a| < 1$, then the opposite is true. For example, if you look at the time series for $a = 0.9$ in figure 7.2, you will see that, after an initial burn in period where the series is affected by the initial condition $x_0$, the process settles down to random motion within a band (between about 5 and 15 in this case).

# Generates an AR(1) time series starting from x = init
ar1ts <- function(init, n, b, a) {
    x <- numeric(n)
    x[1] <- init 
    w <- rnorm(n-1)
    for (t in 1:(n-1)) {
        x[t+1] <- b + a * x[t] + w[t]
    }   
    return(x)
}   

avec <- c(0.1, -0.1, 0.9, -0.9, 1.1, -1.1)
N <- 200
par(mfrow=c(3,2))  # Arrangement of figures 
for (a in avec) {
    plot(ar1ts(0, N, 1, a), type="l", xlab=paste("a = ", a), ylab="")
}

If the scalar Gaussian AR(1) model in (7.3) is unable to capture important aspects of a time series that we wish to model, we can generalize it in several different directions. For example, we can remove the assumption that the shocks are normally distributed, in which case it is simply called the scalar AR(1) model. This might be desirable, if, say, the state process has heavy tails.

We can also generalize the scalar AR(1) model to $\RR^K$, yielding the vector AR(1) model, or VAR(1) model

Here $\{ \boldw_t \}$ is assumed to be IID and to satisfy $\EE \boldw_t = \boldzero$ and $\EE[\boldw_t\boldw_t'] = \boldI$. If $\boldw_t$ is multivariate normal, then (7.4) is called the Gaussian VAR(1). The vector $\boldx_t$ is called the state vector. A huge variety of models can be represented as a VAR(1) process or contain VAR processes as subcomponents (e.g., linear state space models).

An alternative generalization of the scalar AR(1) model is the scalar AR($p$) model, where the next state $x_{t+1}$ is a linear function not just of the current state $x_t$, but of the last $p$ previous states. For example, the AR(2) process has dynamics

Although $x_{t+1}$ is a function of two lagged states, $x_t$ and $x_{t-1}$, we can reformulate it as a first-order model. To begin, we define an additional state variable $y_t$ via $y_t = x_{t-1}$. The dynamics can then be expressed as

We can write this in matrix form as

This is a special case of the VAR(1) model in (7.4). This example shows how higher order processes can be reduced to first-order processes by increasing the number of state variables. For this reason most theoretical discussion below takes place in the context of first-order models.

7.1.2.1Nonlinear Models¶

Linear models are simple but that isn’t always an advantage. Sometimes they are unable to capture the kinds of dynamics we observe in data. For example, the bursts of volatility in the graph of daily returns on the Hang Seng Index from 1990 to 2014 shown in figure 7.3 are difficult to replicate without introducing nonlinearities into our model. It’s also true that many theoretical modeling exercises produce models that are nonlinear.

One well-known nonlinear model is the $p$th order autoregressive conditional heteroskedasticity model (ARCH($p$) model). The model arose from the need to attach a time varying volatility component $\sigma_t$ to asset return models. Suppose, for example, that $x_t = \sigma_t w_t$ where $\sigma_{t+1}^2 = \alpha_0 + \alpha_1 x_t^2$. Combining these equations gives the ARCH(1) model

and $\alpha_0 > 0$, $\alpha_1 \geq 0$. Still better fits to asset returns data can be obtained via Generalized ARCH models, the simplest of which is the GARCH(1,1) process

Now next period volatility depends on its own lagged state as well as $x_t$.

Another popular nonlinear model is the smooth transition autoregression (STAR) model

where $g$ is of the form

Here $\tau \colon \RR \to [0,1]$ is increasing and satisfies $\lim_{s \to -\infty} \tau(s)=0$ and $\lim_{s \to \infty}\tau(s) =1$. For $s$ near $-\infty$, $g(s)$ is similar to $b_0 + a_0 s$. For $s$ near $+\infty$, $g(s)$ is close to $b_1 + a_1 s$. Thus the dynamics transition between two different linear models, with the smoothness of the transition depending on the shape of $\tau$.

7.2Markov Processes¶

All of the examples in Section 7.1.2 are special cases of a general class of process called Markov processes. Markov processes are very rewarding to study because the concepts are broad enough to include a vast array of different stochastic processes, while at the same time having enough structure to support elegant and powerful theory. Here we cover the foundations, from definitions and basic properties to stability and asymptotics.

7.2.1The Markov Assumption¶

Let $S$ be a subset of $\RR^K$ and let $\bB(S)$ be the Borel subsets of $S$. The fundamental primitive of a discrete time Markov process on $S$ is a stochastic kernel or transition probability function, which is a function $Q \colon S \times \bB(S) \to [0, 1]$ such that

$S$ is called the state space of the model. Part (ii) is a mild regularity condition that ensures certain integrals are well defined. Part (i) is the major defining characteristic, yielding a distribution supported on $S$ for every point $\bolds$ in $S$. We think of $Q$ and some initial distribution $P_0$ as generating a corresponding Markov process $\{\boldx_t\}$ taking values in $S$ via

• draw $\boldx_0$ from $P_0$

• for $t$ in $0, 1, 2, \ldots$ do

  • draw $\boldx_{t+1}$ from the distribution $Q(\boldx_t, \cdot)$

• end

The value $Q(\bolds, B)$ is then the probability that $\boldx_{t+1}$ lands in $B$ given $\boldx_t = \bolds$. The sequence $\{\boldx_t\}$ so generated is also called a sample path for $Q$.

There is a measure-theoretic procedure for building this process $\{\boldx_t\}$ more formally but the construction is usually unnecessary because most Markov processes are introduced via a stochastic difference equation, which specifies both the process $\{\boldx_t\}$ and, implicitly, the associated stochastic kernel $Q$. The canonical stochastic difference equation for a first-order Markov process takes the form

Here $\{\boldw_t\}_{t \geq 1}$ is an IID sequence of $\RR^M$-valued shocks with common distribution $\Psi$ and $G$ is a given $\bB$-measurable function mapping the current state $\boldx_t \in S$ and shock $\boldw_{t+1} \in \RR^M$ into the new state $\boldx_{t+1} \in S$. The initial condition $\boldx_0$ and the shocks $\{\boldw_t\}_{t \geq 1}$ are also assumed to be independent of each other. The corresponding stochastic kernel is

The choice of date $t+1$ is irrelevant here because of the IID assumption, and is used only to tie in with (7.7). Recalling that $\Psi$ is the distribution of $\boldw_{t+1}$, the stochastic kernel can also be written as

This expression makes it clear that $Q$ is defined by the two primitives $G$ and $\Psi$.

Returning to the general process (7.7), repeated substitution gives

and so on. Continuing in this fashion, we see that the state vector $\boldx_t$ can be written as a function of $\boldx_0$ and the shocks $\boldw_1,\ldots,\boldw_t$ for any $t$. In other words, for each $t$, there exists a function $H_t$ such that

Often there is no neat expression for the function $H_t$. Nonetheless, (7.8) clarifies the fact that (7.7) pins down each $\boldx_t$ as a well-defined random vector, depending on the initial condition and the shocks up until date $t$. Since $\boldx_t$ is a function of $\boldx_0, \boldw_1, \boldw_2,\ldots, \boldw_t$, and these random vectors are all, by the IID assumption, independent of $\boldw_{t+k}$ for any $k \geq 1$, fact 5.1.10 tells us that the current state $\boldx_t$ and future shocks are independent. Let’s record this:

Although it’s rarely made explicit, behind all the random vectors $\{\boldx_t\}$ lies a single probability space $(\Omega, \fF, \PP)$. The idea is that an element $\omega$ of $\Omega$ is selected by “nature” at the start of time. This determines the initial condition $\boldx_0$ and the shocks $\boldw_t$ as

From these, each state vector $\boldx_t$ is determined via

where $H_t$ is the function in (7.8).

7.2.1.1The Density Case¶

If $Q(\bolds, \cdot)$ is absolutely continuous for all $\bolds$, then we can define the corresponding stochastic density kernel or the transition density via

In other words, $q(\bolds, \cdot)$ is the conditional density of $\boldx_{t+1}$ given $\boldx_t = \bolds$. Heuristically, $q(\bolds, \bolds') \diff \bolds'$ is the probability of transitioning from $\bolds$ to $\bolds'$ in one step.

A visualization of the ARCH kernel (7.11) is given in figure 7.4 when $\alpha_0 = 0.5$ and $\alpha_1 = 0.8$. The diagonal line is the 45 degree line. Each vertical fiber of the graph corresponds to one density $s' \mapsto q(s, s')$. For large $s$ (to the right on the horizontal axis), the corresponding density puts most mass below the 45 degree line, indicating a tendency to decrease. Conversely, when $-s$ is large, the process tends to increase.

7.2.2Marginal and Joint Distributions¶

Let’s note some features of the marginal and joint distributions of a Markov process. Consider our generic process $\boldx_{t+1} = G(\boldx_t, \boldw_{t+1})$ from (7.7), where $\lL(\boldx_0) = P_0$ and $\{\boldw_t\}$ is IID on $\RR^M$ with common distribution $\Psi$. The marginal distribution of the state vector is just $\lL(\boldx_t)$. Since $\boldx_t$ is a well-defined random vector (see in particular (7.8)), $\lL(\boldx_t)$ is also well-defined, and we denote it in what follows by $P_t$. That is,

Figure 7.5 shows a visualization of this sequence when $\mu_0=0$, $\sigma_0 = b = c =1$ and $a = 0.8$. In fact we don’t actually need to recompute the sums in the definition of $\mu_t$ and $\sigma_t^2$ at each step because these sequences satisfy the recursions

You can verify this by taking expectations and variances on both side of (7.3).

For a general Markov process, tracking evolution of the first two moments is usually more complicated, and even then fails to provide full information on the distributions. We do, however, have the following general result:

Exercise 7.6.4 asks you to prove this. Fact 7.2.2 is a version of the law of total probability . To calculate the probability of reaching $B$ next period, we sum the probability of that event given outcome $\bolds$ today, weighted by the probability that $\bolds$ occurs today. Here’s a parallel result for the density case.

Let’s now look at joint distributions. We’ll cover these only in the density case, which is the most important case in econometrics.

This result can be obtained in a few different ways. One is to repeatedly apply (5.27) to factor the joint density into

This is the expression for the joint density of a general stochastic process, where the current state depends on lagged states of all orders. With the first-order Markov property, an expression like $p(\bolds_{t+1} \given \bolds_t, \ldots, \bolds_0)$ reduces to $p(\bolds_{t+1} \given \bolds_t)$, which is just $q(\bolds_t, \bolds_{t+1})$.

7.2.3Stationarity of Markov Processes¶

Consider the notion of stationarity introduced in Section 7.1.1. Not every Markov process is stationary. As stated above, a necessary condition is that the marginals are constant, and in figure 7.5 we saw an example of how the marginals $\{P_t\}$ can change over time. Figure 7.6 shows the same sequence but now longer, with identical parameters. Probability mass moves from left to right.

While $\{P_t\}$ is not typically constant, for the particular sequence shown in figure 7.6, the difference between successive distributions is diminishing over time. We can see this analytically as well by recalling that $P_t = \nN(\mu_t, \sigma^2_t)$, where $\mu_t$ and $\sigma^2_t$ are as given in (7.14). Provided that $|a| < 1$, these sequences converge, with

Let $P_\infty$ denote the distribution associated with these limits. That is,

What happens if we actually attain this limiting distribution, or if we take our initial draw $x_0$ from this distribution? On an intuitive level, it seems that the distribution sequence will be constant, remaining at $P_\infty$ over time. Indeed this is the case, as exercise 7.6.1 asks you to confirm.

The limiting distribution is an example of what is called a stationary distribution. More generally, a stationary distribution for a Markov process on $S$ with stochastic kernel $Q$ is any distribution $P_\infty$ on $S$ satisfying

Comparing with the recursion $P_{t+1}(B) = \int Q(\bolds, B) P_t(\diff \bolds)$ from fact 7.2.2, we see that the definition of stationarity of $P_\infty$ means precisely that

From this we see that any process starting at $P_\infty$ will be identically distributed. In fact it is stationary, in the sense that (7.1) holds. Let’s record this:

Moments of the stationary distribution are called the stationary moments of the Markov process. Other names for the stationary distribution include the "invariant" or "unconditional" distribution, and the moments are then called the invariant or unconditional moments.

Such a density is called a stationary density.

7.2.3.1Existence¶

Not every Markov process has a stationary distribution. For example, in the linear scalar AR(1) model $x_{t+1} = b + a x_t + c w_{t+1}$, we know that the marginal variance evolves according to $\sigma^2_{t+1} = a^2 \sigma_t^2 + c^2$. If $|a| \geq 1$, then this sequence diverges. Because the variance is not constant, the distribution sequence cannot be constant either.

Let’s state one sufficient condition for existence. The condition we give slightly specializes a famous result known as the Krylov–Bogolyubov theorem. It uses the notion of a coercive function, which is any nonnegative function $V$ on the state space $S$ such that

is a compact subset of $S$ for all $\gamma \in \RR$. Intuitively, coercive functions are functions that go to $+\infty$ as we move "away from the center" of the state space.

Now we can state the following:

See §11.2 of Stachurski (2009) for proofs and references. Condition (7.20) is an example of a drift condition. One way to understand it is to write it as

As we move "away from the center" of the state space, $V(\bolds)$ becomes large, and the term on the right is eventually $< 1$. Looking at the left hand side, this tells us that the value attached to the state by $V$ is expected to decrease, which in turn suggests movement away from the "edges" of the state space and back towards the center.

7.2.4Asymptotics of Markov Processes¶

When do Markov processes exhibit nice properties such as ergodicity? To gain some intuition, let’s first consider the process $x_{t+1} = a x_t + w_{t+1}$, where $\{w_t\}$ is IID $\nN(0, 1)$ and $|a| < 1$. We consider a stationary process, where each $x_t$ has marginal

A law of large numbers for this process takes the form

Fact 6.1.2 tells us that for (7.23) to hold it suffices that

The first part is obvious, since $\EE \bar x_T = \frac{1}{T} \sum_{t=1}^T \EE x_t$ and $\EE x_t = 0$. For the second part, since $\EE \bar x_T = 0$,

Using $\EE[x_n x_m] = \EE[x_m x_n]$ we obtain the bound

Thus, to obtain $\var \bar x_T \to 0$ as $T \to \infty$, we require the covariance terms $\EE[x_{t+j} x_t]$ to converge to zero in $j$ sufficiently quickly. That they do in this particular instance is the subject of exercise 7.6.3.

To summarize, if correlations die out sufficiently quickly, then we can regain the LLN. The same is true for the CLT when second moment conditions hold. We now state a much more general result for sample paths of Markov processes. The proofs are omitted but the essential ideas are the same.

Conditions (a) and (b) are not dissimilar to the conditions for existence in theorem 7.2.1. Condition (b) is the drift condition in (7.20) specialized to the density setting. The major difference is that $q$ is required to be strictly positive. This gives enough mixing to generate uniqueness and the sample path properties. Theorem 7.2.2 can be proved from results in Meyn & Tweedie (2009).^[1]

Theorem 7.2.2 is but one of a vast range of results on asymptotic properties of Markov processes. Discussion of alternatives is given in Section 7.5. In fact the conditions are rather strict and can certainly be weakened at the expense of a more complex theorem statement. Nonetheless, they suffice for many applications of interest.

The constant $\sigma^2_h$ in (7.27) is given by the expression

with the understanding that $\{\boldx_t\}$ is a stationary process. If $\sigma^2_h = 0$, then the limiting distribution is degenerate.

To see why, recall from (7.13) that the transition density satisfies $q(\bolds, \bolds' ) = \psi(\bolds' - g(\bolds))$. Since $\psi$ is assumed continuous and positive on $\RR^K$, the transition density has these properties on $\RR^K \times \RR^K$. Hence condition (a) in theorem 7.2.2 is satisfied. Condition (b) is also satisfied because, taking $\lambda$ and $L$ as in (7.29), we have

For most Markov processes the stationary density has no known closed-form. However, in the case of example 7.2.11, the stationary density is known to have the form $p_{\infty}(s) = 2 \phi(s) \Phi(q s)$, where $q := a (1 - a^2)^{-1/2}$, $\phi$ is the standard normal density and $\Phi$ is the standard normal CDF.

7.2.5The Linear Case¶

Recall the VAR process in (7.4), which, to repeat, takes the form

where $\{ \boldw_t \}$ is assumed to be IID and to satisfy $\EE \boldw_t = \boldzero$ and $\EE[\boldw_t\boldw_t'] = \boldI$. Regarding the asymptotics of this process, if the shock vector has a continuous and positive density, then we can apply fact 7.2.7 to obtain asymptotic stability and LLN and CLT results whenever the matrix norm of $\boldA$ satisfies $\|\boldA \| < 1$. To see this, let $g(\bolds) = \boldb + \boldA \bolds$. By the norm triangle inequality and the definition of the matrix norm, we have

Hence (7.29) holds.

As mentioned above, the stability and sample path conditions considered above are relatively strict. In the linear case, we can certainly do better. Here’s one such result:

Here $\rho(\boldA)$ is the spectral radius of $\boldA$, as defined. Under this condition, a unique stationary distribution exists. Taking expectations across (7.30), we see that its mean $\boldmu_{\infty}$ must satisfy $\boldmu_{\infty} = \boldb + \boldA \boldmu_{\infty}$. Applying the Neumann series lemma, this equation has the unique solution

Let $\boldSigma_\infty$ be the asymptotic variance (i.e., the variance–covariance matrix of the stationary distribution). From (7.30) it must satisfy

This is an example of a discrete time Lyapunov equation. Under the condition $\rho(\boldA) < 1$ it can be solved by iteration or by using existing numerical solvers. Here’s an example in Julia that uses a solver from the QuantEcon package:

using QuantEcon

A = [0.8 -0.2; 
    -0.1 0.7]
C = [0.5 0.4;
     0.4 0.6]

solve_discrete_lyapunov(A, C * C')

7.3Martingales¶

Martingales generalize the notion of a random walk, where increments to the current state are draws from a white noise process. The concept of martingales has been central to many advances in probability since the mid twentieth century. Econometricians also like martingales because they neatly capture the idea of "unpredictable" change.

7.3.1Definitions¶

To define martingales, we need the notion of a filtration, which is an increasing sequence of information sets. Recall from Section 5.2.4 that an information set is just a set of random variables. Let $\{\fF_t\}$ be a sequence of information sets. That is, $\fF_t$ is an information set for each $t$. The sequence $\{\fF_t\}$ is called a filtration if, in addition, $\fF_t \subset \fF_{t+1}$ for all $t$. Intuitively, expanding information sets reflect the idea that more information is revealed over time.

Now let $\{m_t\}$ be a scalar stochastic process and let $\{\fF_t\}$ be a filtration. We say that $\{m_t\}$ is adapted to $\{\fF_t\}$ if $m_t$ is $\fF_t$-measurable for every $t$. (For the definition of measurability, see Section 5.2.4.)

The intuitive meaning is this: $\fF_t$ represents the variables we know as of time $t$. If $\{m_t\}$ is adapted to $\fF_t$, then we can compute $m_t$ at time $t$ as well.

The proof goes like this: By adaptedness, $m_t$ is $\fF_t$-measurable. From the definition of a filtration we have $\fF_t \subset \fF_{t+j}$. Applying fact 5.2.10, it follows that $m_t$ is $\fF_{t+j}$-measurable. Fact 7.3.1 now follows from fact 5.2.13.

A sequence of integrable random variables $\{m_t\}$ adapted to a filtration $\{\fF_t\}$ is called a martingale with respect to $\{\fF_t\}$ if

The expression $\EE [ m_{t+1} \given \fF_t]$ is sometimes simplified to $\EE_t [ m_{t+1}]$ when the filtration is understood.

7.3.1.1Martingale Difference Sequences¶

A stochastic process $\{d_t\}$ is called a martingale difference sequence (or MDS) with respect to some filtration $\{\fF_t\}$ if

An MDS is so named because if $\{m_t\}$ is a martingale with respect to $\{\fF_t\}$, then $d_t = m_t - m_{t-1}$ is an MDS with respect to $\{\fF_t\}$. See exercise 7.6.7.

Indeed, by the law of iterated expectations, $\EE d_t = \EE[ \EE[ d_t \given \fF_{t-1} ]] = \EE 0 = 0$.

7.3.2Martingale Difference LLN and CLT¶

Next we consider asymptotics for martingale difference sequences. Such sequences are good candidates for the LLN and CLT. To see this, recall that integrable IID sequences always satisfy the LLN. With finite second moment they also satisfy the CLT. Suppose now that $\{m_t\}$ is a stationary MDS with respect to some filtration $\{\fF_t\}$. We don’t have independence in this case but we do have the following:

To see this, fix $t \geq 0$ and $j \geq 1$. We have

Since $t + j - 1 \geq t$ and $\{\fF_t\}$ is a filtration, we know that $m_t$ is $\fF_{t+j-1}$-measurable. Hence

The lack of correlation in fact 7.3.3 is a key part of the next result.

The proof of (7.32) is identical to that of the classical LLN in theorem 6.1.1. Theorem 7.3.1 can be established from a martingale CLT proved in Durrett (2010).

7.4Simulation¶

Simulation is invaluable as a tool for exploring statistical concepts and solving applied problems. For example, the current popularity of Bayesian methods for econometrics and machine learning has been driven to a large extent by clever simulation algorithms that have proved successful at solving hard estimation problems.

Modern scientific programming environments all have good routines for simulating independent observations from a range of common named distributions. We’ve seen examples of this in the displayed computer code. But it’s easy to run into situations where the distribution we are working with is nonstandard, or has no closed-form solution. We then need to roll our own routines for producing random variates. Let’s review two of the most important algorithms.

7.4.1Inverse Transforms¶

Suppose that we wish to draw from a distribution on $\RR$ with CDF $F$. Assuming that we can generate uniform random variables, we can do so via the inverse transform method:

• draw $u$ from the uniform distribution on $[0, 1]$

• return $F^{-1}(u)$, where $F^{-1}$ is the quantile function of $F$ (see (4.35))

We wish to show that this technique is valid, or, equivalently, that

In the proof we’ll assume that $F$ is strictly increasing. (This isn’t necessary, but the proof of the general case provides no additional intuition.) Let $\lL(u) = U[0, 1]$ and let $z := F^{-1}(u)$. To see that $\lL(z) = F$, observe that, by strict monotonicity, $F(a) \leq F(b)$ if and only if $a \leq b$. As a result, for any $s \in \RR$,

The defining property of a $U[0, 1]$ random variable $u$ is that $\PP\{u \leq t\} = t$ whenever $0 \leq t \leq 1$. Using this fact in the preceding calculation gives $\PP\{z \leq s\} = F(s)$, proving (7.34).

Here’s one extension for the basic inverse transform method. Suppose that we want to simulate from a joint distribution on $\RR^N$ that is defined in terms of a copula (see Section 5.1.5). In particular, suppose that $F_n$ is a continuous CDF on $\RR$ for $n=1,\ldots, N$, that $C$ is a copula on $\RR^N$ and that the joint distribution of interest is

(This is just (5.22).) If we can simulate from $C$, then we can simulate from $F$. The method is

• draw $u_1, \ldots, u_N$ from $C$

• return $(x_1, \ldots, x_N) := (F^{-1}_1 (u_1), \ldots, F^{-1}_N(u_N))$

Exercise 7.6.13 asks you to prove that $(x_1, \ldots, x_N)$ is a draw from $F$.

7.4.2Markov Chain Monte Carlo¶

Markov chain Monte Carlo (MCMC) is a way of simulating from a given density $\pi$ on $S \subset \RR^K$. The idea is to construct a stochastic kernel $P$ on $S$ such that

See (7.26) for a clearer statement of (ii). If we specialize to indicator functions we can approximate probabilities under $\pi$ too.

There are two major strands of MCMC algorithms: Metropolis–Hastings and the Gibbs sampler. Here we’ll cover the former. While boldface symbols are omitted in what follows, the ideas are valid in any dimension. In the section on Bayesian estimation in Section 8.3.4, I will explain why Metropolis–Hastings fits so beautifully with Bayesian updating.

The Metropolis–Hastings algorithm starts with a Markov process in the form of a stochastic density kernel $q = q(s, s')$ called the proposal density. Draws from the proposal density are called proposals. Each time we draw a proposal, we either (i) accept it and move to that new state, or (ii) reject it and stay put. The probability of accepting is structured so that the chain tends to stay in regions where $\pi$ puts most probability mass. As a consequence, for a sequence $\{x_t\}$ generated by this process, we get

This is a version of (7.26) with $h = \1_B$.

The ideas are subtle, so let’s build up in steps. Consider first a general accept-reject style process of the kind discussed above, where the probability of acceptance is determined by an arbitrary function $\alpha = \alpha(x_t, y)$ of the current state $x_t$ and proposal $y$. The function takes values in $[0, 1]$. The algorithm for generating $x_{t+1}$ from $x_t$ given $\alpha$ and proposal density $q$ is given in algorithm 1.

The stochastic kernel $P$ associated with $\{x_t\}$ in algorithm 1 has the form

for all $s \in S$ and $B \in \bB(S)$, where

To show this, we can break the probabilities up as follows. Fix any Borel set $B$ and let $x_t$ be given at $s$. Let $x_{t+1}$ be drawn according to algorithm 1 when $x_t = s$. By additivity of $\PP$ we can write $\PP\{x_{t+1} \in B\}$ as

Our aim is to show that the two terms in (7.36) match the two terms on the right-hand side of (7.35).

Using the fact that intersections become products when we switch to indicators (fact 4.1.4, page ), we can applying the law of iterated expectations to obtain

The term $\EE \left[ \1 \{u \leq \alpha(s, y)\} \given y \right]$ is just $\PP\{u \leq \alpha(s, y)\}$ with $y$ given, which, since $u$ is uniform, evaluates to $\alpha(s, y)$. Hence the final expression in the last display evaluates to

Equation (7.37) says that the first term in (7.36) is equal to the first term on the right-hand side of (7.35). For the second terms, note that $x_{t+1} = x_t = s$ when $u > \alpha(s, y)$, so

To match this expression with the second term in (7.35), we just need to show that $\PP \{u \leq \alpha(s, y)\} = \int p(s, s') \diff s'$. But this is immediate from (7.37) if we take $B = S$.

So far we have verified that the kernel $P$ in (7.35) is the stochastic kernel corresponding to algorithm 1. Now let’s double back and design the selection function $\alpha$ so that the stationary distribution of $P$ is the target density $\pi$. In the Metropolis–Hastings algorithm, this is done by setting