If you're going to approximate an AR(1) process using a finite state Markov-chain, many people will use the method given in either Tauchen (1986) or Tauchen and Hussey (1991). For a highly persistent process, however, Rouwenhorst (1995) outperforms both of these methods considerably. A rigorous comparison is given by Kopecky and Suen (2010).

Python code for Tauchen and Hussey is in a previous post. An implementation of Rowenhorst is given here (and shown below).

Python code for Tauchen and Hussey is in a previous post. An implementation of Rowenhorst is given here (and shown below).

import numpy as np
def rouwen(rho, mu, step, num):
'''
Adapted from Lu Zhang and Karen Kopecky. Python by Ben Tengelsen.
Construct transition probability matrix for discretizing an AR(1)
process. This procedure is from Rouwenhorst (1995), which works
well for very persistent processes.
INPUTS:
rho - persistence (close to one)
mu - mean and the middle point of the discrete state space
step - step size of the even-spaced grid
num - number of grid points on the discretized process
OUTPUT:
dscSp - discrete state space (num by 1 vector)
transP - transition probability matrix over the grid
'''
# discrete state space
dscSp = np.linspace(mu -(num-1)/2*step, mu +(num-1)/2*step, num).T
# transition probability matrix
q = p = (rho + 1)/2.
transP = np.array([[p**2, p*(1-q), (1-q)**2], \
[2*p*(1-p), p*q+(1-p)*(1-q), 2*q*(1-q)], \
[(1-p)**2, (1-p)*q, q**2]]).T
while transP.shape[0] <= num - 1:
# see Rouwenhorst 1995
len_P = transP.shape[0]
transP = p*np.vstack((np.hstack((transP,np.zeros((len_P,1)))), np.zeros((1, len_P+1)))) \
+ (1-p)*np.vstack((np.hstack((np.zeros((len_P, 1)), transP)), np.zeros((1, len_P+1)))) \
+ (1-q)*np.vstack((np.zeros((1, len_P+1)), np.hstack((transP, np.zeros((len_P, 1)))))) \
+ q*np.vstack((np.zeros((1, len_P+1)), np.hstack((np.zeros((len_P, 1)), transP))))
transP[1:-1] /= 2.
# ensure columns sum to 1
if np.max(np.abs(np.sum(transP, axis=1) - np.ones(transP.shape))) >= 1e-12:
print('Problem in rouwen routine!')
return None
else:
return transP.T, dscSp

Hey guys,

ReplyDeleteJust a comment to thank you for your blog, it's almost a daily reading to me. Keep going, we are several researchers at my lab to love it!

Best,

Maxime

Thanks Maxime! We'll be sure to keep posting!

ReplyDeleteShould an exception be raised if "if num % 2 == 1"?

ReplyDeleteHi Brian, thanks for pointing that out. Turns out it isn't necessary.

Deletewhy do u take a transpose at the end?

ReplyDeleteJust a convention I prefer.

Delete