Computational Economics

Go to top, next, previous, or johnstachurski.net

Stochastic Recursive Sequences

Previously we studied finite Markov chains

Next we consider stochastic dynamic models on infinite state spaces

The presentation is terse: See EDTC, chapter 6 for details

The Model

Let's focus on the class of models given by:

where

Each X_t takes values in a set S
- X_t is called the state variable, and S is the state space
Each W_t takes values in a set Z
- W_t is called the shock, and Z is the shock space
ψ is a distribution on S
- called the initial condition
φ is a distribution on Z
- called the shock distribution

These kinds of models are called

stochastic recursive sequences (SRSs), or
stochastic difference equations, or
(first order, time homogenous) Markov chains (or Markov processes)

Example: Stochastic Solow--Swan Growth Model

The law of motion for the capital stock is given by

where

with state and shock spaces S = Z = (0,∞)

Note that in this model, the savings rate is fixed

In modern economic analysis, savings is usually chosen optimally via dynamic programming

We'll get to that soon

Example: STAR model

The smooth transition threshold autoregression (STAR) model is of the form

where

The state and shock spaces S and Z are just the real line, and

is a smooth transition function such as a Gaussian cumulative distribution function

Both of these models are nonlinear.

More difficult to analyze than linear systems
- But more interesting too
Analysis is aided by simulation and other computational techniques

Simulation and Time Series

Listing 6.1 in EDTC provides a class for simulating time series of an abstract SRS

## Filename: srs.py
## Author: John Stachurski

class SRS:
 
    def __init__(self, F=None, phi=None, X=None):
        """Represents the model X_{t+1} = F(X_t, W_{t+1}), where W ~ phi

        Parameters: 
            * F and phi are functions, where phi() returns a draw from phi 
            * X is a number representing the initial condition
        """
        self.F, self.phi, self.X = F, phi, X

    def update(self):
        "Update the state according to X = F(X, W)."
        self.X = self.F(self.X, self.phi())

    def sample_path(self, n):
        "Generate path of length n from current state."
        path = []
        for i in range(n):
            path.append(self.X)
            self.update()
        return path

Problem:

Use this code to generate 5 time series (sample paths) for the Solow model above when

For the parameters use

alpha, sigma, s, delta = 0.5, 0.2, 0.5, 0.1

The initial conditions for k₀ should range over a grid of length 5 from 1 to 80 (one for each time series)

Each series should be of length 500

Plot the resulting series on a single graph as follows

Solution:

Here's the Solow model

## Filename: solmod.py
## Author: John Stachurski

import numpy as np

alpha, sigma, s, delta = 0.5, 0.2, 0.5, 0.1  

def F(k, z): 
    return s * (k**alpha) * z + (1 - delta) * k

def phi():
    return np.random.lognormal(sigma=sigma)

And here's the solution to the exercise

## Filename: solowts.py
## Author: John Stachurski

import srs 
import solmod
import numpy as np
import matplotlib.pyplot as plt

solow_srs = srs.SRS(F=solmod.F, phi=solmod.phi)
initial_conditions = np.linspace(1, 80, num=5)
for k in initial_conditions:
    solow_srs.X = k
    plt.plot(solow_srs.sample_path(500))

plt.show()

We've used two files in order to separate the description of the data (i.e., model) from the logic of the solution

Marginal Distributions

Consider again the general SRS

Suppose we are interested in what happens at some fixed point in time T

The random variable X_T is defined recursively by the model

We can generate independent observations of X_T as follows

Draw an observation of X₀ from its distribution ψ
Iterate forward up until time T
Record the value of X_T
Repeat (i.e., go to step 1)

If we loop n times then we get n independent observations of X_T

These are draws from the marginal distribution of X_T

Problem:

Add a method marginal to the SRS class which implements this algorithm

Parameters are init, T and n
Returns a NumPy array of n observations of X_T when X₀ = init

Solution:

## Filename: srs2.py
## Author: John Stachurski

import numpy as np

class SRS:
    
    def __init__(self, F=None, phi=None, X=None):
        """ Represents the model X_{t+1} = F(X_t, W_{t+1}), where W ~ phi

        Parameters: 
            * F and phi are functions, where phi() returns a draw from phi 
            * X is a number representing the initial condition
        """
        self.F, self.phi, self.X = F, phi, X

    def update(self):
        "Update the state according to X = F(X, W)."
        self.X = self.F(self.X, self.phi())

    def sample_path(self, n):
        "Generate path of length n from current state."
        path = np.zeros(n)
        for t in range(n):
            path[t] = self.X
            self.update()
        return path

    def marginal(self, init=None, T=None, n=None):
        "Returns n draws of X_T, starting from the state X=init."
        samples = np.zeros(n)
        for i in range(n):
            self.X = init
            for t in range(T):
                self.update()
            samples[i] = self.X
        return samples

Here's an example of usage

import srs2
import solmod   # Contains the definition of the Solow model (in solmod.py)
import pylab as pl
solow_srs = srs2.SRS(F=solmod.F, phi=solmod.phi)
ts = solow_srs.marginal(init=1, T=10, n=5000)
pl.hist(ts, bins=30, normed=True, label="s = 0.5")
solmod.s = 0.6  # Reset the savings rate from 0.5 (default) to 0.6
ts = solow_srs.marginal(init=1, T=10, n=5000)
pl.hist(ts, bins=30, normed=True, label="s = 0.6")
solmod.s = 0.7  # Reset the savings rate from 0.6 to 0.7
ts = solow_srs.marginal(init=1, T=10, n=5000)
pl.hist(ts, bins=30, normed=True, label="s = 0.7")
pl.legend()
pl.show()

The plot looks like this

Can you give some interpretation of this figure?

Nonparametric Kernel Density Estimates

As we saw, the time T distribution can be represented as a histogram

If the density to be estimated is smooth, there are better ways

In particular, we can use the nonparametric kernel density estimator (NPKDE)

The general expression of the NPKDE (for the density f of a real-valued random variable X) is

In this expression

K is some density (often, the standard normal density)
and h is either a parameter or a function of the data
- usually referred to as the bandwidth

The estimate is a collection of n "bumps," one centered on each data point

These are then summed and normalized to create a density

The bandwidth parameter plays a role similar to the number of bins used in the histogram:

high value means that the densities we place on each data point are flat with large tails
A low value means they are concentrated around each data point, and estimate is spiky

Common default for the bandwidth is to use Silverman's rule

where σ is the standard deviation of the sample

Silverman's rule is optimal when the density to be estimated is normal (i.e., Gaussian)

Problem:

Implement the nonparametric kernel density estimator

Write as a class, where data, bandwidth and kernel are the instance variables
Supply the standard normal pdf as the default kernel
Use Silverman's rule to give a default bandwidth
Provide a __call__ method which evaluates the estimate at a point x

Solution:

## Filename: npkde.py
## Author: John Stachurski

import numpy as np

def g(x):
    return (1.0 / np.sqrt(2 * np.pi)) * np.exp(- 0.5 * x**2)

class NPKDE:

    def __init__(self, X=None, h=None, K=g):
        """
        Nonparametric kernel density estimator.

        Parameters:
            * X is a NumPy array containing the observations
            * h is a number
            * K is a vectorized callable which represents a density 
        """
        self.X, self.h, self.K = X, h, K 
        if self.h == None:
            # If bandwidth is not defined, use Silverman's rule of thumb
            self.h = 1.06 * np.std(X) * len(X)**(-1.0 / 5)

    def __call__(self, x):
        return (1.0 / self.h) * np.mean(self.K((x - self.X) / self.h))

Here's an example of usage

>>> import numpy as np
>>> import npkde
>>> import matplotlib.pyplot as plt
>>> f = npkde.NPKDE(X=np.random.randn(250))
>>> xgrid = np.linspace(-4, 4, num=200)
>>> plt.plot(xgrid, [f(x) for x in xgrid])
>>> plt.show()

The plot looks like this

Problem:

Using the NPKDE, investigate the density of k_T for different savings rates

savings_rates = np.linspace(0.5, 0.9, num=5)

Otherwise the model parameters are as above

The marginal distributions are specified as follows

k0 = 1      # Initial condition
T = 10      # Date of the marginal density to be computed
n = 500     # Number of observations at each date

Plot the marginal densities of k_T for each savings rate on the same graph

Your output should look something like this (add the legend if you can)

Solution:

## Filename: solow_npkde.py
## Author: John Stachurski

import srs2 
import solmod
import npkde
import numpy as np
import matplotlib.pyplot as plt

k0 = 1      # Initial condition
T = 10      # Date
n = 500     # Number of observations at each date

savings_rates = np.linspace(0.5, 0.9, num=5)
solow_srs = srs2.SRS(F=solmod.F, phi=solmod.phi)

# Get all the observations for each savings rate
obs_list = []
for s in savings_rates:
    solmod.s = s
    obs_list.append(solow_srs.marginal(k0, T, n))

# Find the smallest and largest of all observations
kmin = np.array(obs_list).min()
kmax = np.array(obs_list).max()
xgrid = np.linspace(kmin, kmax, num=100)

for i, s in enumerate(savings_rates):
    f = npkde.NPKDE(obs_list[i])
    fig_label = 's = ' + str(s)
    plt.plot(xgrid, [f(x) for x in xgrid], label=fig_label)

plt.legend(loc='best')
plt.show()

Look-Ahead Estimators

Actually there is a better way to compute these marginal distributions

Uses the look-ahead estimator

We saw this estimator earlier, in the case of finite Markov models

Can be applied to estimating marginal and stationary densities of continuous state models such as the stochastic Solow model

See EDTC, section 6.1 for an introduction

The programming is not hard, and we wont discuss it further

Lectures on Computational Economics

John Stachurski

Stochastic Recursive Sequences

The Model

Example: Stochastic Solow--Swan Growth Model

Example: STAR model

Simulation and Time Series

Marginal Distributions

Nonparametric Kernel Density Estimates

Look-Ahead Estimators