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
- X
- Each W
t takes values in a set Z- W
t is called the shock, and Z is the shock space
- W
- ψ 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 k0 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 XT is defined recursively by the model
We can generate independent observations of XT as follows
- Draw an observation of X0 from its distribution ψ
- Iterate forward up until time T
- Record the value of XT
- Repeat (i.e., go to step 1)
If we loop n times then we get n independent observations of XT
These are draws from the marginal distribution of XT
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 XT when X0 = 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 kT 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 kT 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
0 comments:
Post a Comment