One of the first steps in the Matched-Expectations Path (MEP) algorithm for solving Heterogeneous Agent models with Aggregate Shocks is to generate a T-period sequence of values of the Markov process on the aggregate shock. The obvious and simplest way to do this is to simulate the process —pick an arbitrary value for period 1, then use a random number generator together with the markov transition matrix to iteratively create periods 2 to T— which is formally known as Monte Carlo simulation of the markov chain. This is the approach used by both Hanbaek Lee’s codes, and the current implementation of MEP in VFI Toolkit.
Creating this sequence is low cost (relative to the rest of the algorithm) and the choice of sequence is likely important to the performance of the MEP algorithm (as an extreme, there is a non-zero probability we just simulate a path that takes the same value of the aggregate shock in every period in which case the MEP will fail). This suggests that if we can find a way to ‘reliably’ produce the sequence of aggregate shocks with ‘good properties’, and if the cost of doing so is moderate, then it will improve the performance of MEP. Hopefully it will allow using smaller T, but at a minimum it should make the algorithm more ‘stable’ across runs.
In what follows I will first explain MC, QMC, and RQMC in the context of uniformly distributed i.i.d random variables, and then once we have established these concepts I will extend them to more general i.i.d. random variables and to markov chains. To fix ideas, we will initially pretend that the ‘good property’ we want our shock sequences to have is simply an unbiased low-variance estimator of the mean, and then later we will think more generally what ‘good properties’ might mean in the context of the MEP algorithm.
Imagine that our interest is around simulating an i.i.d random variable X \sim Uniform(0,1), and the ‘good property’ we care about is estimating the mean of the random variable, \mu=E[X]. Focusing on a uniform i.i.d. variable will make it easier to introduce the concepts, and will turn out to easily generalize to cover lot’s of other situations.
Monte Carlo (MC) method estimates \mu by
\bar{X}_n = \frac{1}{n} \sum_{i=1}^{n} X_i
where X_1, X_2, \dots, X_n are n independent U(0,1) random variables. To implement this we simply use a random number generator to create the X_i and then compute their sum.
Using random sequences like this leaves us open to the aforementioned issue that on rare occasion we will be very unlucky in our random draw and get, e.g., that all n=10 draws of our U(0,1) are less than 0.1.
Quasi-Monte Carlo (QMC) replaces the independent random draws X_i with a set of deterministic points P_n=\{y_1,y_2,\dots,y_n\} that cover [0,1) more evenly. QMC estimates \mu by,
\bar{\mu}_n = \frac{1}{n} \sum_{i=1}^{n} y_i
Roughly speaking, P_n is called a highly-uniform point set or low-discrepency point set if some measure of discrepency between the empirical distribution of P_n and the uniform distribution converges to 0 faster than \mathcal{O}(n^{-1/2}), which is the typical rate for independent random points, as n \to \infty.
There are two main problems with QMC, both stemming directly from the fact it is deterministic. The first is that QMC will have a non-zero error in the estimate (we cannot call it biased, due to the non-random nature you cannot even define what unbiased would mean but it is a problem in a similar vein to bias). The second is that QMC does not allow us to compute the variance of the estimator (again, inescapably due to the non-random nature of QMC). Both of these are ‘fixed’ by Randomized Quasi-Monte Carlo (RQMC), which as it turns out also happily happens to provide convergence of the estimator at least as fast as QMC, and in some situations even faster than QMC, as well as having a few other minor nicer properties than QMC.
Randomized Quasi-Monte Carlo (RQMC) takes the QMC and shifts/scrambles the points such that
(i) each individual point is uniformly distributed (0,1), marginally
(ii) the set of points collectively retain their low-discrepancy structure
Implementing RQMC involves applying some form of shifting/scrambling that constructs our random set \{z_1,z_2,\dots,z_n\} from the QMC set \{y_1,y_2,\dots,y_n\}, and this shifting/scrambling is done in such a way as to introduce (i), while maintaining (ii). RQMC estimates \mu by,
\bar{\mu}_{n,rqmc} = \frac{1}{n} \sum_{i=1}^{n} z_i
a later post will discuss how to perform shifting/scrambling, but some of the main methods are lattice rules and digital nets. As long as we are careful about both how the deterministic sequence was created, and about how the shifting/scrambling was implemented the RQMC gives us an unbiased and low-variance estimator for \mu. Simulation studies show that RQMC often substantially outperforms both MC and QMC, with the gains depending on the application.
So far everything has been done with a one-dimensional Uniform(0,1) i.i.d. random variable. All of this extends trivially to s-dimensions, with a Uniform(0,1)^s i.i.d random variables. Since most i.i.d. distributions can be expressed as a function of a uniform variable, we simply switch all the above sums to look like \frac{1}{n} \sum_{i=1}^{n} f(x_i) with f(x_i) in place of x_i and this covers essentially all other random variables of interest. [In fact ‘create uniform and transform to desired distribution’ is how computers create draws from essentially all the other distributions.] Notice that this same f(x_i) setting also covers when our interest is in the behaviour of higher-order moments, rather than the mean.
That will do for this post. We have seen the ideas of MC, QMC and RQMC for i.i.d. random variables. In the next post I will describe how to construct the deterministic set for QMC and how to perform the shifting/scrambling involved in getting from QMC to RQMC. The post after that will discuss how to extend all of this to Markov Chains. The post after that will discuss the application to MEP.
PS. This post leans heavily on L’Ecuyer - Randomized Quasi-Monte Carlo: An Introduction for Practitioners.
PPS. Of course in practice the random numbers used by computers are not truly random. Instead they come from a pseudo-random number generator, a deterministic sequence constructed to look random (Matlab’s default is the Mersenne Twister). Good pseudo-random number generators are themselves designed on discrepancy-like criteria: the set of all points the generator can ever produce is made as uniform as possible (1,2). So your “random” numbers are, in a sense, a random slice of one enormous QMC-like point set. There is no theorem that MC with a pseudo-random number generator converges exactly as MC with true randomness would, but modern generators pass extensive empirical test batteries and this has not been found to matter in practice.