quant-econ Solutions: Continuous State Markov Chains

Solutions for http://quant-econ.net/stationary_densities.html



In [1]:

    
%matplotlib inline



In [2]:

    
import numpy as np
import matplotlib.pyplot as plt

Exercise 1

Look ahead estimation of a TAR stationary density, where the TAR model is

$$ X_{t+1} = \theta |X_t| + (1 - \theta^2)^{1/2} \xi_{t+1} $$

and $\xi_t \sim N(0,1)$. Try running at n = 10, 100, 1000, 10000 to get an idea of the speed of convergence.



In [3]:

    
from scipy.stats import norm, gaussian_kde
from quantecon import LAE

phi = norm()
n = 500
theta = 0.8
# == Frequently used constants == #
d = np.sqrt(1 - theta**2) 
delta = theta / d

def psi_star(y):
    "True stationary density of the TAR Model"
    return 2 * norm.pdf(y) * norm.cdf(delta * y) 

def p(x, y):
    "Stochastic kernel for the TAR model."
    return phi.pdf((y - theta * np.abs(x)) / d) / d

Z = phi.rvs(n)
X = np.empty(n)
for t in range(n-1):
    X[t+1] = theta * np.abs(X[t]) + d * Z[t]
psi_est = LAE(p, X)
k_est = gaussian_kde(X)

fig, ax = plt.subplots(figsize=(10,7))
ys = np.linspace(-3, 3, 200)
ax.plot(ys, psi_star(ys), 'b-', lw=2, alpha=0.6, label='true')
ax.plot(ys, psi_est(ys), 'g-', lw=2, alpha=0.6, label='look ahead estimate')
ax.plot(ys, k_est(ys), 'k-', lw=2, alpha=0.6, label='kernel based estimate')
ax.legend(loc='upper left')









    Out[3]:





<matplotlib.legend.Legend at 0x3fca450>

Exercise 2

Here's one program that does the job



In [4]:

    
from scipy.stats import lognorm, beta

# == Define parameters == #
s = 0.2
delta = 0.1
a_sigma = 0.4       # A = exp(B) where B ~ N(0, a_sigma)
alpha = 0.4         # f(k) = k**alpha

phi = lognorm(a_sigma) 

def p(x, y):
    "Stochastic kernel, vectorized in x.  Both x and y must be positive."
    d = s * x**alpha
    return phi.pdf((y - (1 - delta) * x) / d) / d

n = 1000     # Number of observations at each date t
T = 40       # Compute density of k_t at 1,...,T

fig, axes = plt.subplots(2, 2, figsize=(10, 10))
axes = axes.flatten()
xmax = 6.5

for i in range(4):
    ax = axes[i] 
    ax.set_xlim(0, xmax)
    psi_0 = beta(5, 5, scale=0.5, loc=i*2)  # Initial distribution

    # == Generate matrix s.t. t-th column is n observations of k_t == #
    k = np.empty((n, T))
    A = phi.rvs((n, T))
    k[:, 0] = psi_0.rvs(n)
    for t in range(T-1):
        k[:, t+1] = s * A[:,t] * k[:, t]**alpha + (1 - delta) * k[:, t]

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

    ygrid = np.linspace(0.01, xmax, 150)
    greys = [str(g) for g in np.linspace(0.0, 0.8, T)]
    greys.reverse()
    for psi, g in zip(laes, greys):
        ax.plot(ygrid, psi(ygrid), color=g, lw=2, alpha=0.6)
    ax.set_xlabel('capital')



In [ ]: