2. Continuous State Markov Chains#

2.2.1. Definitions and Basic Properties#

In our lecture on finite Markov chains, we studied discrete-time Markov chains that evolve on a finite state space $S$.

In this setting, the dynamics of the model are described by a stochastic matrix — a nonnegative square matrix $P=P[i,j]$ such that each row $P[i,\cdot ]$ sums to one.

The interpretation of $P$ is that $P[i,j]$ represents the probability of transitioning from state $i$ to state $j$ in one unit of time.

In symbols,

Equivalently,

• $P$ can be thought of as a family of distributions $P[i,\cdot ]$, one for each $i\in S$

• $P[i,\cdot ]$ is the distribution of $X_{t+1}$ given $X_t=i$

(As you probably recall, when using NumPy arrays, $P[i,\cdot ]$ is expressed as P[i,:])

In this section, we’ll allow $S$ to be a subset of $\mathbb{R}$, such as

• $\mathbb{R}$ itself

• the positive reals $(0,\mathrm{\infty })$

• a bounded interval $(a,b)$

The family of discrete distributions $P[i,\cdot ]$ will be replaced by a family of densities $p(x,\cdot )$, one for each $x\in S$.

Analogous to the finite state case, $p(x,\cdot )$ is to be understood as the distribution (density) of $X_{t+1}$ given $X_t=x$.

More formally, a stochastic kernel on $S$ is a function $p:S\times S\to \mathbb{R}$ with the property that

1. $p(x,y)\ge 0$ for all $x,y\in S$

2. $\int p(x,y)dy=1$ for all $x\in S$

(Integrals are over the whole space unless otherwise specified)

For example, let $S=\mathbb{R}$ and consider the particular stochastic kernel $p_w$ defined by

What kind of model does $p_w$ represent?

The answer is, the (normally distributed) random walk

To see this, let’s find the stochastic kernel $p$ corresponding to (2.2).

Recall that $p(x,\cdot )$ represents the distribution of $X_{t+1}$ given $X_t=x$.

Letting $X_t=x$ in (2.2) and considering the distribution of $X_{t+1}$, we see that $p(x,\cdot )=N(x,1)$.

In other words, $p$ is exactly $p_w$, as defined in (2.1).

2.2.2. Connection to Stochastic Difference Equations#

In the previous section, we made the connection between stochastic difference equation (2.2) and stochastic kernel (2.1).

In economics and time-series analysis we meet stochastic difference equations of all different shapes and sizes.

It will be useful for us if we have some systematic methods for converting stochastic difference equations into stochastic kernels.

To this end, consider the generic (scalar) stochastic difference equation given by

Here we assume that

• $\{{\xi }_t\}\stackrel{\text{ IID }}{\sim }\varphi$, where $\varphi$ is a given density on $\mathbb{R}$

• $\mu$ and $\sigma$ are given functions on $S$, with $\sigma (x)>0$ for all $x$

Example 1: The random walk (2.2) is a special case of (2.3), with $\mu (x)=x$ and $\sigma (x)=1$.

Example 2: Consider the ARCH model

Alternatively, we can write the model as

This is a special case of (2.3) with $\mu (x)=\alpha x$ and $\sigma (x)=(\beta +\gamma x^2{)}^{1/2}$.

Example 3: With stochastic production and a constant savings rate, the one-sector neoclassical growth model leads to a law of motion for capital per worker such as

Here

• $s$ is the rate of savings

• $A_{t+1}$ is a production shock

  • The $t+1$ subscript indicates that $A_{t+1}$ is not visible at time $t$

• $\delta$ is a depreciation rate

• $f:{\mathbb{R}}_{+}\to {\mathbb{R}}_{+}$ is a production function satisfying $f(k)>0$ whenever $k>0$

(The fixed savings rate can be rationalized as the optimal policy for a particular set of technologies and preferences (see [Ljungqvist and Sargent, 2018], section 3.1.2), although we omit the details here).

Equation (2.5) is a special case of (2.3) with $\mu (x)=(1-\delta )x$ and $\sigma (x)=sf(x)$.

Now let’s obtain the stochastic kernel corresponding to the generic model (2.3).

To find it, note first that if $U$ is a random variable with density $f_U$, and $V=a+bU$ for some constants $a,b$ with $b>0$, then the density of $V$ is given by

(The proof is below. For a multidimensional version see EDTC, theorem 8.1.3).

Taking (2.6) as given for the moment, we can obtain the stochastic kernel $p$ for (2.3) by recalling that $p(x,\cdot )$ is the conditional density of $X_{t+1}$ given $X_t=x$.

In the present case, this is equivalent to stating that $p(x,\cdot )$ is the density of $Y:=\mu (x)+\sigma (x){\xi }_{t+1}$ when ${\xi }_{t+1}\sim \varphi$.

Hence, by (2.6),

For example, the growth model in (2.5) has stochastic kernel

where $\varphi$ is the density of $A_{t+1}$.

(Regarding the state space $S$ for this model, a natural choice is $(0,\mathrm{\infty })$ — in which case $\sigma (x)=sf(x)$ is strictly positive for all $s$ as required)

2.2.3. Distribution Dynamics#

In this section of our lecture on finite Markov chains, we asked the following question: If

1. $\{X_t\}$ is a Markov chain with stochastic matrix $P$

2. the distribution of $X_t$ is known to be ${\psi }_t$

then what is the distribution of $X_{t+1}$?

Letting ${\psi }_{t+1}$ denote the distribution of $X_{t+1}$, the answer we gave was that

This intuitive equality states that the probability of being at $j$ tomorrow is the probability of visiting $i$ today and then going on to $j$, summed over all possible $i$.

In the density case, we just replace the sum with an integral and probability mass functions with densities, yielding

It is convenient to think of this updating process in terms of an operator.

(An operator is just a function, but the term is usually reserved for a function that sends functions into functions)

Let $\mathcal{D}$ be the set of all densities on $S$, and let $P$ be the operator from $\mathcal{D}$ to itself that takes density $\psi$ and sends it into new density $\psi P$, where the latter is defined by

This operator is usually called the Markov operator corresponding to $p$

With this notation, we can write (2.9) more succinctly as ${\psi }_{t+1}(y)=({\psi }_{t}P)(y)$ for all $y$, or, dropping the $y$ and letting “$=$” indicate equality of functions,

Equation (2.11) tells us that if we specify a distribution for ${\psi }_0$, then the entire sequence of future distributions can be obtained by iterating with $P$.

It’s interesting to note that (2.11) is a deterministic difference equation.

Thus, by converting a stochastic difference equation such as (2.3) into a stochastic kernel $p$ and hence an operator $P$, we convert a stochastic difference equation into a deterministic one (albeit in a much higher dimensional space).

2.2.4. Computation#

To learn about the dynamics of a given process, it’s useful to compute and study the sequences of densities generated by the model.

One way to do this is to try to implement the iteration described by (2.10) and (2.11) using numerical integration.

However, to produce $\psi P$ from $\psi$ via (2.10), you would need to integrate at every $y$, and there is a continuum of such $y$.

Another possibility is to discretize the model, but this introduces errors of unknown size.

A nicer alternative in the present setting is to combine simulation with an elegant estimator called the look-ahead estimator.

Let’s go over the ideas with reference to the growth model discussed above, the dynamics of which we repeat here for convenience:

Our aim is to compute the sequence $\{{\psi }_t\}$ associated with this model and fixed initial condition ${\psi }_0$.

To approximate ${\psi }_t$ by simulation, recall that, by definition, ${\psi }_t$ is the density of $k_t$ given $k_0\sim {\psi }_0$.

If we wish to generate observations of this random variable, all we need to do is

1. draw $k_0$ from the specified initial condition ${\psi }_0$

2. draw the shocks $A_1,\dots ,A_t$ from their specified density $\varphi$

3. compute $k_t$ iteratively via (2.12)

If we repeat this $n$ times, we get $n$ independent observations $k_t^1,\dots ,k_t^n$.

With these draws in hand, the next step is to generate some kind of representation of their distribution ${\psi }_t$.

A naive approach would be to use a histogram, or perhaps a smoothed histogram using SciPy’s gaussian_kde function.

However, in the present setting, there is a much better way to do this, based on the look-ahead estimator.

With this estimator, to construct an estimate of ${\psi }_t$, we actually generate $n$ observations of $k_{t-1}$, rather than $k_t$.

Now we take these $n$ observations $k_{t-1}^1,\dots ,k_{t-1}^n$ and form the estimate

where $p$ is the growth model stochastic kernel in (2.8).

What is the justification for this slightly surprising estimator?

The idea is that, by the strong law of large numbers,

with probability one as $n\to \mathrm{\infty }$.

Here the first equality is by the definition of ${\psi }_{t-1}$, and the second is by (2.9).

We have just shown that our estimator ${\psi }_t^n(y)$ in (2.13) converges almost surely to ${\psi }_t(y)$, which is just what we want to compute.

In fact, much stronger convergence results are true (see, for example, this paper).

2.2.5. Implementation#

A class called LAE for estimating densities by this technique can be found in lae.py.

Given our use of the __call__ method, an instance of LAE acts as a callable object, which is essentially a function that can store its own data (see this discussion).

This function returns the right-hand side of (2.13) using

• the data and stochastic kernel that it stores as its instance data

• the value $y$ as its argument

The function is vectorized, in the sense that if psi is such an instance and y is an array, then the call psi(y) acts elementwise.

(This is the reason that we reshaped X and y inside the class — to make vectorization work)

Because the implementation is fully vectorized, it is about as efficient as it would be in C or Fortran.

2.2.6. Example#

The following code is an example of usage for the stochastic growth model described above

# == Define parameters == #
s = 0.2
δ = 0.1
a_σ = 0.4                    # A = exp(B) where B ~ N(0, a_σ)
α = 0.4                      # We set f(k) = k**α
ψ_0 = beta(5, 5, scale=0.5)  # Initial distribution
ϕ = lognorm(a_σ)


def p(x, y):
    """
    Stochastic kernel for the growth model with Cobb-Douglas production.
    Both x and y must be strictly positive.
    """
    d = s * x**α
    return ϕ.pdf((y - (1 - δ) * x) / d) / d

n = 10000    # Number of observations at each date t
T = 30       # Compute density of k_t at 1,...,T+1

# == Generate matrix s.t. t-th column is n observations of k_t == #
k = np.empty((n, T))
A = ϕ.rvs((n, T))
k[:, 0] = ψ_0.rvs(n)  # Draw first column from initial distribution
for t in range(T-1):
    k[:, t+1] = s * A[:, t] * k[:, t]**α + (1 - δ) * k[:, t]

# == Generate T instances of LAE using this data, one for each date t == #
laes = [LAE(p, k[:, t]) for t in range(T)]

# == Plot == #
fig, ax = plt.subplots()
ygrid = np.linspace(0.01, 4.0, 200)
greys = [str(g) for g in np.linspace(0.0, 0.8, T)]
greys.reverse()
for ψ, g in zip(laes, greys):
    ax.plot(ygrid, ψ(ygrid), color=g, lw=2, alpha=0.6)
ax.set_xlabel('capital')
ax.set_title(f'Density of $k_1$ (lighter) to $k_T$ (darker) for $T={T}$')
plt.show()

The figure shows part of the density sequence $\{{\psi }_t\}$, with each density computed via the look-ahead estimator.

Notice that the sequence of densities shown in the figure seems to be converging — more on this in just a moment.

Another quick comment is that each of these distributions could be interpreted as a cross-sectional distribution (recall this discussion).

2.3.1. Example and Definitions#

To illustrate the issues, recall that Hopenhayn and Rogerson [Hopenhayn and Rogerson, 1993] study a model of firm dynamics where individual firm productivity follows the exogenous process

As is, this fits into the density case we treated above.

However, the authors wanted this process to take values in $[0,1]$, so they added boundaries at the endpoints 0 and 1.

One way to write this is

If you think about it, you will see that for any given $x\in [0,1]$, the conditional distribution of $X_{t+1}$ given $X_t=x$ puts positive probability mass on 0 and 1.

Hence it cannot be represented as a density.

What we can do instead is use cumulative distribution functions (cdfs).

To this end, set

This family of cdfs $G(x,\cdot )$ plays a role analogous to the stochastic kernel in the density case.

The distribution dynamics in (2.9) are then replaced by

Here $F_t$ and $F_{t+1}$ are cdfs representing the distribution of the current state and next period state.

The intuition behind (2.14) is essentially the same as for (2.9).

2.3.2. Computation#

If you wish to compute these cdfs, you cannot use the look-ahead estimator as before.

Indeed, you should not use any density estimator, since the objects you are estimating/computing are not densities.

One good option is simulation as before, combined with the empirical distribution function.

2.4.1. Theoretical Results#

Analogous to the finite case, given a stochastic kernel $p$ and corresponding Markov operator as defined in (2.10), a density ${\psi }^{\ast }$ on $S$ is called stationary for $P$ if it is a fixed point of the operator $P$.

In other words,

As with the finite case, if ${\psi }^{\ast }$ is stationary for $P$, and the distribution of $X_0$ is ${\psi }^{\ast }$, then, in view of (2.11), $X_t$ will have this same distribution for all $t$.

Hence ${\psi }^{\ast }$ is the stochastic equivalent of a steady state.

In the finite case, we learned that at least one stationary distribution exists, although there may be many.

When the state space is infinite, the situation is more complicated.

Even existence can fail very easily.

For example, the random walk model has no stationary density (see, e.g., EDTC, p. 210).

However, there are well-known conditions under which a stationary density ${\psi }^{\ast }$ exists.

With additional conditions, we can also get a unique stationary density ($\psi \in \mathcal{D}\text{ and }\psi =\psi P⟹\psi ={\psi }^{\ast }$), and also global convergence in the sense that

This combination of existence, uniqueness and global convergence in the sense of (2.16) is often referred to as global stability.

Under very similar conditions, we get ergodicity, which means that

for any (measurable) function $h:S\to \mathbb{R}$ such that the right-hand side is finite.

Note that the convergence in (2.17) does not depend on the distribution (or value) of $X_0$.

This is actually very important for simulation — it means we can learn about ${\psi }^{\ast }$ (i.e., approximate the right-hand side of (2.17) via the left-hand side) without requiring any special knowledge about what to do with $X_0$.

So what are these conditions we require to get global stability and ergodicity?

In essence, it must be the case that

1. Probability mass does not drift off to the “edges” of the state space.

2. Sufficient “mixing” obtains.

For one such set of conditions see theorem 8.2.14 of EDTC.

In addition

• [Stokey et al., 1989] contains a classic (but slightly outdated) treatment of these topics.

• From the mathematical literature, [Lasota and MacKey, 1994] and [Meyn and Tweedie, 2009] give outstanding in-depth treatments.

• Section 8.1.2 of EDTC provides detailed intuition, and section 8.3 gives additional references.

• EDTC, section 11.3.4 provides a specific treatment for the growth model we considered in this lecture.

2.4.2. An Example of Stability#

As stated above, the growth model treated here is stable under mild conditions on the primitives.

• See EDTC, section 11.3.4 for more details.

We can see this stability in action — in particular, the convergence in (2.16) — by simulating the path of densities from various initial conditions.

Here is such a figure.

All sequences are converging towards the same limit, regardless of their initial condition.

The details regarding initial conditions and so on are given in this exercise, where you are asked to replicate the figure.

2.4.3. Computing Stationary Densities#

In the preceding figure, each sequence of densities is converging towards the unique stationary density ${\psi }^{\ast }$.

Even from this figure, we can get a fair idea what ${\psi }^{\ast }$ looks like, and where its mass is located.

However, there is a much more direct way to estimate the stationary density, and it involves only a slight modification of the look-ahead estimator.

Let’s say that we have a model of the form (2.3) that is stable and ergodic.

Let $p$ be the corresponding stochastic kernel, as given in (2.7).

To approximate the stationary density ${\psi }^{\ast }$, we can simply generate a long time-series $X_0,X_1,\dots ,X_n$ and estimate ${\psi }^{\ast }$ via

This is essentially the same as the look-ahead estimator (2.13), except that now the observations we generate are a single time-series, rather than a cross-section.

The justification for (2.18) is that, with probability one as $n\to \mathrm{\infty }$,

where the convergence is by (2.17) and the equality on the right is by (2.15).

The right-hand side is exactly what we want to compute.

On top of this asymptotic result, it turns out that the rate of convergence for the look-ahead estimator is very good.

The first exercise helps illustrate this point.