Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
91,816 characters · 28 sections · 77 citation commands
A Scrambled Method of Moments
JEL Classification: C11, C12, C13, C32, C36.\newline Keywords: Simulated Method of Moments, Indirect Inference, quasi-Monte Carlo, Scramble.
\baselineskip=18.0pt \thispagestyle{empty} \setcounter{page}{0}
Simulation-based estimation is a popular approach to estimate complex economic models. The econometrician simply matches sample with simulated moments, drawn from a model of interest. The resulting Simulated Method of Moments (SMM) or Indirect Inference estimator makes estimation feasible even though the likelihood or the moments' expectation, required for MLE and GMM, may be impossible or impractical to compute. However, using simulations rather than analytical computations introduces simulation noise, which increases the variance of the estimates. Also, the resulting simulated objective function is typically non-smooth and hence more difficult to minimize numerically. In theory, simulating many samples reduces simulation noise and smoothes the objective function, making optimization easier. In practice, however, this may not be feasible because of the increased computational cost. Also, more informative moments used for Indirect Inference or the Efficient Method of Moments can to be more computationally demanding than simpler moments used in SMM, leading to a computational tradeoff between the informativeness of the moments and simulation noise. A practical solution is to use variance reduction techniques such as antithetic draws\footnote{See Section (ref) for a brief overview of antithetic sampling.} which can reduce simulation noise with nearly no computational overhead.
This paper investigates quasi-Monte Carlo (qMC) integration, another variance reduction approach, in the context of simulation-based estimation, with an emphasis on the scramble of owen1995randomly. Under certain conditions, qMC can approximate an expectation at a faster rate than the usual Monte-Carlo (MC) Central Limit Theorem (CLT). This suggests that a Scrambled Method of Moments could outperform conventional SMM estimates using as many simulated samples. This is shown to be the case for a large class of models in cross-sections and short-panels with potentially non-smooth moments as in McFadden1989,Pakes1989 or auxiliary parameters as in Gourieroux1993. For time series, qMC may not apply directly because of a curse of dimensionality over the time dimension. A class of models and moments which circumvent this issue are described.
Using the scramble in an estimation setting poses several practical and theoretical challenges. These sequences are designed to approximate a fixed integral of an iid sequence. Improper use of the scramble under dependence or with covariates may result in inconsistent estimators. Hence, the first and main contribution of the paper is methodological. The second contribution is theoretical. Uniform Laws of Large Numbers (ULLN) and CLTs are provided to handle smooth moments in cross-sections and short-panels. Scrambled draws are random and identically distributed but not independent. This makes it more challenging to handle time series and non-smooth moments. A stochastic equicontinuity result for cross-sections and short-panels is established by re-writing the scrambled empirical process as the sum of a non-identically distributed but independent array with a standard qMC sequence. This allows to invoke existing results for each term separately. In the time series setting, a similar strategy allows to invoke results for bounded dependent heterogeneous arrays.
The finite sample properties of the Scrambled Method of Moments are illustrated using several simple Monte-Carlo examples including an income process with \scalebox{-1}[1]{”} lots of heterogeneity” BROWNING2010. In this example, the scramble improves on SMM with random and antithetic draws in terms of variance. Furthermore, optimization over the 2,000 replications was on average 15% faster with the scramble than SMM using as many simulated samples - because scrambled moments are smoother than MC moments.
After a review of the literature, Section (ref) provides an overview of (quasi)-Monte Carlo integration which is lesser known in economics. Section (ref) shows how to implement the Scrambled Method of Moments in various settings. Asymptotic results for each algorithm are given in Section (ref) and the proofs are in Appendix (ref). Section (ref) illustrates its finite sample properties using Monte-Carlo simulations. Section (ref) concludes.
There are two related literatures: simulation-based estimation and variance reduction techniques. In economics, simulation-based estimation includes the Simulated Method of Moments McFadden1989,Pakes1989,duffie-singleton, Indirect Inference Gourieroux1993 and the Efficient Method of Moments Gallant1996. See smith-palgrave for an overview of simulation-based estimation in economics and common empirical applications. In statistics, Bayesian methods such as Approximate Bayesian Computation Marin2012 and Synthetic Likelihood Wood2010 are more common. See jjng-14 for an overview and comparisons of these frequentist and Bayesian methods.
As discussed in the introduction, Monte-Carlo methods introduce simulation noise in the estimation which increases the variance of the estimator. There is a large number of variance reduction techniques, the following summary will only cover some of those that are most relevant to simulation-based estimation. One approach is to use low-discrepancy sequences - this is more commonly known as quasi-Monte Carlo integration. These sequences were initially designed to compute integrals of iid sequences and can achieve faster than $\sqrt{n}$-rate convergence. More details are given in Section (ref). qMC integration has been extended to non-linear state-space filtering Gerber2015,Gerber2017, MCMC sampling Owen2005 and importance sampling for ABC estimation Buchholz2017. A key takeaway from these papers is that a lot of care is required in implementing qMC integration in non iid settings (MCMC or filtering) where `naive' implementations may be inconsistent. This may explain why it is only rarely used in empirical economics, even though their appeal has been known for some time judd1998numerical. In economics, antithetic draws are a popular variance reduction technique. However, they can lead to either efficiency gains or losses depending on the integrand as discussed in Section (ref). Another variance reduction method, which is more popular in statistics, is the control variates approach robert2013monte.\footnote{Note that despite the similarity in names, this is not related to control variable estimation used in structural econometric estimation.} The main idea is to augment the estimating sample and simulated moments with analytically tractable moments for the shocks themselves. This additional information can help reduce the uncertainty attributable to simulation noise.\footnote{See Davis2019 for an application of control variates to Indirect Inference. Control variates were also considered for qMC integration in Hickernell2005a.} The control functional approach Oates2017, which uses all the information about the distribution of the shocks, can result in faster than $\sqrt{n}$-rate convergence. Important efficiency gains require the control variate moments to be sufficiently rich which could lead to a curse of dimensionality. For instance, the model of Section (ref) has shocks with dimension $d=30$ so that spanning polynomials of order up to $2$ or $3$ would require introducing $496$ or $5,456$ additional moments respectively. The number of moments quickly becomes greater than the sample size itself.
The following provides a brief overview of Monte-Carlo (MC) and quasi-Monte Carlo (qMC) integration.\footnote{For further reading, Lemieux2009 provides a non-technical introduction to MC and qMC integration; Dick2010 provide the underlying theory.} Throughout, we are interested in evaluating the integral of a known measurable function $f: [0,1]^d \to \mathbb{R}$:
by using a fixed or random sequence of points $u_1,\dots,u_n$ in $[0,1]^d$:
A widely applicable approach is MC integration. Take iid uniform draws $u_i \sim \mathcal{U}_{[0,1]^d}$ and compute the sample analog $\hat I^{MC}_n = \frac{1}{n}\sum_{i=1}^n f(u_i).$ Assuming $f(u_i)$ has finite variance, $\hat I^{MC}_n$ is unbiased and the approximation error $|\hat I^{MC}_n-I|$ is of order $\sqrt{\text{var}[f(u_i)]/n}$. This implies that in order to reduce the approximation error tenfold, the number of draws must be a hundred times greater: the computational cost increases faster than the approximation error declines.
A popular variance reduction approach is to use antithetic draws. For $n$ even, compute: \[ \hat I_n^{Anti} = \frac{1}{n} \sum_{i=1}^{n/2}[f(u_i) + f(1-u_i)], u_i \overset{iid}{\sim} \mathcal{U}_{[0,1]^d}.\] This approach is only valid if $f(u_i)$ and $f(1-u_i)$ have the same distribution; for instance, $e_i = \Phi^{-1}(u_i) \sim \mathcal{N}(0,1)$ and $-e_i = \Phi^{-1}(1-u_i) \sim \mathcal{N}(0,1)$ as well. Without this property, when the distribution is asymmetric, $\hat I_n^{Anti}$ may not be consistent for $I$.
Assuming $f(u_i)$ and $f(1-u_i)$ have the same distribution, $\hat I_n^{Anti}$ is unbiased and $\text{var}(\hat I_n^{Anti}) = \left( \text{var}[f(u_i)] + \text{cov}[f(u_i),f(1-u_i)] \right)/n$. If $\text{corr}[f(u_i),f(1-u_i)] = -1$, then $\text{var}(\hat I_n^{Anti}) =0$; the estimator is exact as soon as $n=2$. This improves significantly on MC integration. However, if $\text{corr}[f(u_i),f(1-u_i)] = +1$ then $\text{var}(\hat I_n^{Anti}) = 2 \text{var}(\hat I_n^{MC})$. Now, $\hat I_n^{MC}$ outperforms $\hat I_n^{Anti}$.
The performance of antithetic draws relative to simple MC draws will typically depend on both the parameter of interest and the choice of estimating moments. To illustrate, consider the following two examples. First, suppose $I=\mathbb{E}(e_i)$ where $e_i = \Phi^{-1}(u_i) \sim \mathcal{N}(0,1)$. Note that $-e_i \sim \mathcal{N}(0,1)$ and $\hat I_n^{Anti}$ is consistent for $I$ in this example. Since $\text{corr}(e_i,-e_i) = -1$, $\text{var}(\hat I_n^{Anti})=0$. The estimator is exact as soon as $n=2$. Second, suppose $I=\mathbb{E}(e_i^2)$ with $e_i$ as above. Now, $\text{corr}( e_i^2,[-e_i]^2 )=+1$ and $\text{var}(\hat I_n^{Anti}) = 2 \text{var}(\hat I_n^{MC})$. These examples suggest that the moments need to have some asymmetry properties in order to produce efficiency gains. This can be hard to check for intractable non-linear models.
The discussion above shows that some sequences can outperform MC integration. For instance, for $f$ smooth and $u_i \in [0,1]^d$ with $d=1$ the lattice sequence $u_i = i/(n-1), i=0,\dots,n-1$, the estimator \[ \hat I_n^{Lattice} = \frac{1}{n} \sum_{i=1}^{n}f(u_i), \quad u_i = i/(n-1),\, i \in \{0,\dots,n-1\}\] has an approximation error of order $O(\|\partial_u f\|_{\infty}/n)$. The approximation error declines linearly with the computational cost. However for $d\geq 2$, this sequence has approximation errors of order $n^{-1/d}$ which is worse than MC as soon as $d \geq 3$.
It is possible to break this curse of dimensionality. To achieve this, the qMC literature relies on two pivotal inequalities. The first one is the Koksma-Hlawka inequality:
where $\|f\|_{TV}$ is the total variance norm of $f$ in the sense of Hardy and Krause:
$\partial^{|\mathfrak{u}|} f (\mathfrak{u})/\partial \mathfrak{u}$ consists of all univariate derivatives $\partial_{u_1} f(u),\dots,\partial_{u_d} f(u)$ and partial cross-derivatives $\partial^2_{u_1,u_2} f(u),\partial^2_{u_1,u_3} f(u),\partial^2_{u_2,u_3} f(u),\dots,\partial^2_{u_{d-1},u_d} f(u)$ up to order $d$ with $\partial^d_{u_1,\dots,u_d} f(u)$. It does not include repeated derivatives such as $\partial^2_{u_1,u_1} f(u)$. What matters here is the smoothness of $f$ across the co-ordinates $u_1,\dots,u_d$. As a result, integrating over larger dimensions $d$ typically requires additional smoothness in $f$ over these cross-derivatives.
The other term in the Koksma-Hlawka inequality is $D_n^\star(u_1,\dots,u_n)$ which corresponds to the star discrepancy of the sequence $(u_1,\dots,u_n)$, defined as:
In statistics, this is known as the Kolmogorov-Smirnov (KS) distance between the empirical CDF of $(u_1,\dots,u_n)$ and the population CDF of a uniform $\mathcal{U}_{[0,1]^d}$ distribution.
For a given function $f$, reducing the approximation error in ((ref)) implies finding sequences with smaller $D_n^\star$. For iid random draws, $D_n^* = O_p(n^{-1/2})$ by Donsker theorems VanderVaart1996. The lattice sequence above has $D_n^* = O(n^{-1/d})$ for $d \geq 1$.
For any sequence $(u_1,\dots,u_n)$, the second pivotal inequality - initially due to Roth1954 and generalized by Schmidt1970 - provides a lower bound on its star discrepancy:
where $C_d$ is a universal constant which only depends on the dimension $d$. Note the striking difference with the Discrepancy of the set $D_n(u_1,\dots,u_n) = \sup_{i = 1,\dots,n}(\inf_{j \neq i}\|u_i-u_j\|)$ which cannot decrease faster than $n^{-1/d}$. Under the sup-norm distance there is the well known curse of dimensionality which affects grid searches, non-parametric estimation, etc. Under the KS distance, this lower bound suggests that the impact of dimensionality is much less severe.\\
For any $d \geq 1$ fixed, this lower bound suggests a faster than $\sqrt{n}$-rate is achievable. There are a number of deterministic sequences which are close in rate to the bound ((ref)); these include the Sobol, Halton, van der Corput and Hammersley sequences. Most of these are readily available in statistical softwares.\footnote{The R package randtoolbox, the Sobol module in Julia, the quasirandomset toolset in Matlab and the SamplePack library in C++ can generate the Sobol sequence, for instance.} There is, however, a caveat that when $d$ becomes large $C_d$ can also become large. For instance, $C_d = 2^d$ for the Sobol sequence which increases very rapidly with $d$. As a result, in finite samples MC integration may outperform qMC integration for $d$ large relative to $n$. Under additional smoothness conditions, so-called higher-order sequences can achieve even faster rates of order $n^{-\alpha}$ (up to log-terms) for some $\alpha>1$ which depends on the smoothness of $f$ and higher-order properties of the sequence.
These results above are restrictive since they require the integrand $f$ to be smooth, otherwise $\|f\|_{TV} = +\infty$ and the Koksma-Hlawka inequality ((ref)) is uninformative. In economics, many problems involve non-smooth integrands such as simulation-based estimation of discrete choice models Train2009. Also, $\hat I_{n}^{qMC}$ computed with a deterministic sequence is typically biased and its approximation error is hard to evaluate numerically. This would make it difficult to compute standard errors in an estimation setting.
One solution is to use randomized quasi-Monte Carlo (RqMC) methods. A simple randomizer is the digital shift. Take one random draw $u \sim \mathcal{U}_{[0,1]^d}$, a qMC sequence $u_1,\dots,u_n$ ( e.g. Sobol, Halton) and compute $\tilde u_i = [u_i + u] \text{ modulo } 1$. The modulo operator is applied one dimension at a time. This shifts all the co-ordinates of $u_1,\dots,u_n$ by the same random quantity $u$ and preserves the order of magnitude of its star discrepancy $D_n^\star$. The randomized $\tilde u_i$ are identically distributed $\mathcal{U}_{[0,1]^d}$ but not independent. The estimator $\hat I_n^{RqMC}$ is unbiased. To approximate its variance, apply the digital shift with different draws $u$ to compute the integral several times and then compute the variance across these estimates Lemieux2009.
Another randomization approach, which will be the main focus of this paper, is the scramble introduced by owen1995randomly. Similarly to the random shift above, it transforms a deterministic low-discrepancy sequence into random identically but not independently distributed uniform $\mathcal{U}_{[0,1]^d}$ draws. Since the scrambled draws are uniform, the estimator $\hat I_n^{scramble}$ is unbiased. The procedure is described in the box below. The scramble does not deteriorate the discrepancy of the original sequence, in fact it was shown that it can further improve it Dick2010.
The scramble approximates $I$ under the same conditions as the classical CLT as shown in Theorem (ref). The underlying theory is quite involved since it relies of Walsh expansions, an approach similar to Fourier expansions but in a digital basis $b$ which requires an understanding of both number theory and functional approximation theory. See Dick2010 for an introduction to the relevant material and proofs.
Under additional smoothness conditions $\hat I_n^{scramble}$ approximates $I$ at a near $n^{-3/2}$-rate; which is faster than deterministic qMC sequences. A refinement of the orginal algorithm, higher-order scrambling, can achieve even faster rate for smooth integrands; for instance, in some cases the convergence can be of order $n^{-5/2}$ or $n^{-7/2}$. The scrambled estimator $\hat I_n^{Scramble}$ is unbiased. Its variance can be approximated the same way as for $\hat I_n^{RqMC}$. Note that these results assume $d \geq 1$ is fixed. In practice, MC may outperform the scramble for $d$ large. Other scrambles have also been proposed by Hickernell1996 and Matousek1998, among others. See Lemieux2009 and Dick2010 for additional references.\\
Figure (ref) illustrates the differences between random, deterministic and scrambled sequences. A particular realization of a random sample may have points clustered in some areas of $[0,1]^2$, as visible in the figure. The Sobol sequence covers the square better in this case, though some points might appear to cluster here as well. The two realizations of the scrambled Sobol sequence cover the square well and cluster slightly less (better discrepancy) than the deterministic Sobol points.
This section introduces the main algorithms to implement the Scrambled Method of Moments and Scrambled Indirect Inference. The data generating process (DGP) is the same as in Gourieroux1993:
A simple transformation allows to replace $u_{i,t}$ with $e_{i,t} = \Phi^{-1}(u_{i,t}) \overset{iid}{\sim} \mathcal{N}(0,1)$ or other distributions by the Rosenblatt transform. $i=1,\dots,n$ indexes individuals and $t=1,\dots,T$ the time dimension. $y_{i,t}$ is the vector of observed outcome variables. $x_{i,t}$ is a vector of strictly exogenous covariates and $z_{i,t}$ a vector of unobserved latent variables. The functions $g_{obs}$ and $g_{latent}$ are assumed to be known up to a finite dimensional parameter $\theta$ to be estimated.
For static models, which correspond to cross-sections and short-panels, the $t$ index will be omitted to re-write ((ref))-((ref)) and the moments as:
where $y_i = (y_{i,1},\dots,y_{i,T})$ and $u_{i} = (u_{i,1},\dots,u_{i,T}) \in [0,1]^{T \times d_{u_{i,t}}}$. The dimension $d$ of $u_i$ is $T \times \text{dim}(u_{i,t})$, using the notation in ((ref))-((ref)). Given a vector of moments $\hat \psi_n$ and a weighting matrix $W_n$, a simple SMM estimator is given in Algorithm (ref).
Without covariates $x_i$, the expectation $\mathbb{E}[\hat \psi_n^S(\theta)]$ has the same form as ((ref)). The scramble can be applied if the moments have finite variance. The resulting Algorithm (ref) is thus very similar to SMM.
In practice, one samples an $(nS) \times d$ matrix of scrambled shocks rather than $S$ different $n \times d$ matrices of random numbers. This is may be useful because using a large simulated sample of $n \times S$ observations implies a reduction in variance greater than $S$, as a consequence of the faster rate in Theorem (ref), compared to using $S$ independent simulated samples. Asymptotic results for Algorithm (ref) are provided in Proposition (ref), assuming the moments are smooth in $\theta$. These assumptions are comparable to those required for SMM.
When $\hat \psi_n$ is a vector of auxiliary moments as Gourieroux1993, the results from Proposition (ref) can be extended for the scramble as shown in Proposition (ref). Again, the assumptions are comparable to those required for Indirect Inference. These Indirect Inference results could also be extended to non-smooth moments and time series given appropriate changes to the assumptions.
In the presence of covariates, $\mathbb{E}[\hat \psi_n^S(\theta)]$ does not have the same form as ((ref)): \[ \mathbb{E}[\hat \psi_n^S(\theta)] = \int_{[0,1]^d \times \mathcal{X}} \psi\left( g(x,u;\theta),x\right) f_x(x)dx du, \] where $f_x$ is joint density of the covariates $x$.\footnote{The results could be extended to allow some components of $x$ to be discrete. However, the assumptions in the next Section imply that at least one of the covariates should have a continuous density.} Without further assumptions, it is typically not possible to sample from the population $f_x$ directly so that qMC sequence with $n \times S$ elements for $(x_i,u_i)$ cannot be constructed. Taking the covariates as given, Algorithm (ref) relies on $S$ independent scrambled sequences of size $n$ rather than a large sequence of size $n \times S$ as in Algorithm (ref).\footnote{It is implicitly assumed that $(x_1,\dots,x_n)$ is a random sample. If the ordering is deterministic, then $x_i$ and $u_i)$ are not independent. Randomly shuffling the covariates without replacement solves this issue.}
For each $s \in \{1,\dots,S\}$, the function $\mathbb{E}[ \psi( y^s_i(\theta),x_i)|\tilde u_i^s = u]$ does not depend on $x$ so that Theorem (ref) can be applied to this conditional expectation, assuming it has finite variance. This insight was used to derive CLTs for moments based on hybrid sequences which combine MC draws with qMC sequences in Okten2006 and Buchholz2017 for bounded $\psi$. The results in Proposition (ref) extend these results to unbounded empirical processes over $\theta \in \Theta$, allowing $\hat \psi_n^s$ to be non-smooth in $\theta$. The assumptions are more demanding than for SMM, although they could be weakened for smooth moments with covariates. The conditional expectation $\mathbb{E}[\hat \psi_n^S(\cdot)|\tilde u_1,\tilde u_2,\dots]$ itself is required to be smooth in $\theta$, i.e. integrating out the covariates smoothes out the sample and simulated moments. This implies that at least one of the covariates has a continuous density.
For dynamic models, which correspond to time series observations, the $i$ index will be omitted to re-write ((ref))-((ref)) and the moments as:
Covariates are omitted to simplify the theoretical results. Only moments involving a fixed and finite number of lags $L$ will be considered as explained below. Algorithm (ref) details the SMM procedure to estimate ((ref))-((ref)).
To understand the issues caused by the dynamics for the scramble and qMC integration, note that for any initial value $(y_0,z_0)$, $y_t$ can be re-written as: \[ y_t = g_t(u_t,\dots,u_1,y_0,z_0;\theta), \] for some function $g_t$ which can be expressed in terms of $g_{obs}$ and $g_{latent}$. Using this notation, the expected value of $\hat \psi_T$ can be re-written as: \[ \mathbb{E}(\hat \psi_T) = \frac{1}{T} \sum_{t=L+1}^T \int_{[0,1]^{t \times d}} \psi \circ (g_t,\dots,g_{t-L})(u_t,\dots,u_1,y_0,z_0) du_t\dots du_1. \] The expectation above differs from the qMC setting in ((ref)) in several ways. First, the function to be integrated involves $g_t$ which varies with $t$ unlike the function in ((ref)). Second, the integral is computed over $u_1,\dots,u_t$ which has a dimension $t$ that increases with the sample size. This implies a curse of dimensionality for qMC which requires the dimension $d$ to be fixed. Third, both randomized and non-randomized qMC sequences are identically but not independently distributed. A naive implementation of the scramble could introduce spurious dependence in the simulated data and the resulting estimator may not be consistent as a result.
Implementing qMC integration in a dynamic setting without additional structure comes at a cost. In finance, qMC sequences are used to simulate long time series and price financial derivatives Paskov1995,Lemieux2009. This is done by setting $d = T$ and sampling a very large number $n$ of financial series. In the present setting, this amounts to picking $S$ very large and $d=T$ which is not computationally attractive compared to standard SMM.\footnote{Recall that for the Sobol sequence $C_d = 2^d$ so that the error would be of the order of $2^T/S$. Consistency of the qMC integral would require $S \gg 2^T$, i.e. $S$ needs to grow exponentially fast with the sample size $T$.} For state-space filtering, Gerber2015,Gerber2017 propose a Hilbert sorting step to re-sample draws into a low-discrepancy sequence using the Hilbert fractal map from $[0,1]$ to $[0,1]^d$. This Hilbert map can be challenging to implement in practice and suffers from a curse of dimensionality.
The class of moments described in ((ref)) where the number of lags $L$ is fixed and finite allows to circumvent these issues. To get some intuition, suppose that it is possible to draw $(y^1_t,z^1_t) = F_{y,z}^{(-1)}(v^1_t)$ from the stationary distribution directly using the Rosenblatt transform with $v^1_t \sim \mathcal{U}_{[0,1]^{\text{dim}(v_t^1)}}$. Then, using additional shocks $u_t^2,\dots,u_t^L$ one could simulate a short time series consisting of $L \geq 1$ observations for each $t = 1,\dots,T \times S$:
The resulting draws $(y_t^1,\dots,y_t^L)$ are iid over $t =1,\dots,T\times S$ from the stationary distribution by construction.\footnote{This idea was also used in Davis2019 but as a variance reduction method with MC draws.} This is now within the setting of ((ref)). Algorithm (ref) describes a Scrambled Method of Moments for models where simulating as described above is feasible. The main idea is to simulate the $(y_t^1,\dots,y_t^L)$ $T \times S$ times with scrambled shocks $(v_t^1,u_t^2,\dots,u_t^L)_{t = 1,\dots,T\times S} \in [0,1]^d$ with dimension $d=\text{dim}(v_t^1,u_t^2,\dots,u_t^L)$. which depends on the dimension of the shocks and the numbers of lags $L$. Note that while Algorithm (ref) requires $T \times S$ draws, Algorithm (ref) effectively requires $n \times S \times L$ draws. However, the latter Algorithm is massively parallel over $t$ so that for some models it may run faster than the former in a parallel environment. Proposition (ref) provides the asymptotic results for Algorithm (ref).\footnote{When simulating the initial draw with the Rosenblatt transform is not possible, one may consider using a fixed starting value and a burn-in period assuming some decay conditions hold. This is only considered for the hybrid MC-qMC method, theoretical investigations for qMC-only draws is left to future research.}
Sampling from the stationary distribution directly is feasible for some DGPs such as the Gaussian ARMA model (see the Monte-Carlo example in Section (ref)) or the following stochastic volatility process: \[ \log(\sigma_t) = \mu_\sigma + \rho_\sigma \log(\sigma_{t-1}) + \kappa_\sigma e_{t,1}, \quad y_t = \sigma_t e_{t,2},\quad (e_{t,1},e_{t,2}) \overset{iid}{\sim} \mathcal{N}(0,I_2). \] Since the log-volatility follows a Gaussian AR(1) process, one can simply draw $\log(\sigma_t^1) \sim \mathcal{N}(\mu_\sigma/(1-\rho_{\sigma},\kappa_\sigma^2/[1-\rho_\sigma^2])$ and $y_t^1 = \sigma^1_t e^1_{t,2}$ where $e^1_{t,2}$ to simulate $(y_t^1,\sigma_t^1)$ from their stationary distribution. For more complex DGPs this may not be feasible, however.
When the direct approach in Algorithm (ref) is not feasible, an alternative is to sample the initial draws $(y_t^1,z_t^1)$ by MC methods and then simulate $(y_t^2,z_t^2),\dots,(y_t^L,z_t^L)$ using the scramble. This hybrid MC-qMC approach allows to sample from intractable distributions while retaining some of the the features of qMC integration.
The resulting Algorithm (ref) combines elements from Algorithms (ref) and (ref). It requires an additional loop compared to the latter, which is more computationally demanding. Because the estimation combines MC with qMC, the variance of the estimates will typically be greater than a qMC only approach. Note that once the $(z_t^1(\theta),y_t^1(\theta))$ are drawn by MC simulations, $(z_t^\ell(\theta),y_t^\ell(\theta))_{\ell > 1}$ can be simulated in parallel which can be computationally attractive.
Proposition (ref) provides asymptotic results for Algorithm (ref) with conditions similar to duffie-singleton but assuming bounded moments. Relaxing this assumption would require to extend existing CLTs for dependent heterogeneous arrays White1984 which goes beyond the scope of this paper. The simulations in Section (ref) suggest that the estimator performs well with unbounded moments in practice.
Given that the scramble is different from standard Monte-Carlo methods, the following shows how to compute standard errors for $\hat \theta_n^{S}$ for SMM, antithetic draws and the scramble.
Under regularity conditions, the Simulated and Scrambled Method of Moments estimators satisfy the following asymptotic expansion: \[ \hat \theta_n^{S} - \theta_0 = - \left( G^\prime W_n G \right)^{-1} G^\prime W_n \left[\hat \psi_n - \hat \psi_n^S(\theta_0) \right] + o_p(n^{-1/2}), \] where $G = \partial_\theta \mathbb{E} \left[ \hat \psi_n^S(\theta_0) \right]$ is the usual Jacobian matrix. Under a CLT, the asymptotic variance is given by the usual sandwich formula. Given that $W_n$ is chosen by the user, only two terms need to be approximated: the Jacobian $G$ and the asymptotic variance of $[\hat \psi_n - \hat \psi_n^S(\theta_0)]$.
When the moments are smooth, the plug-in Jacobian estimator $\hat G_n = \partial_\theta \hat \psi_n^S(\hat \theta_n^S)$ is consistent for $G$ under a ULLN. For non-smooth moments, there are several possibilities. The more computationally demanding approach is to Bootstrap the estimator $\hat \theta_n^S$ directly. Alternatively, Bruins2018 propose to smooth the draws $y_{i,t}^s$ in dynamic discrete choice models using a kernel; this transforms non-smooth and unbiased into smooth but biased simulated moments. Frazier2019 rely on a change of variable argument to compute analytical Jacobians in a class of discrete choice models. The quasi-Jacobian matrix in Forneron2019 smoothes the moments themselves to approximate $G$. It is also possible to use MCMC methods to sample from a quasi-posterior distributions which approximates the frequentist distribution of $\hat \theta_n^S$ chernozhukov-hong,Wood2010.
For cross-sections and short panels, the asymptotic variance of $[\hat \psi_n - \hat \psi_n^S(\theta_0)]$ in SMM can be approximated with the cross-sectional variance of $[\psi(y_i,x_i)-\frac{1}{S} \sum_{s=1}^S \psi(y_i^s(\hat \theta_n^S),x_i)]$. Pooling all the simulated samples that way ensures that the estimator is consistent for both standard and antithetic draws.\footnote{Another approach is to use use the variance of $\psi(y_i^s(\hat \theta_n^S),x_i)$ divided by $S$ as an estimate for $\hat \psi_n^S(\theta_0)$. Although commonly used, this may actually not be consistent in the presence of antithetic draws. Depending on the correlation described in Section (ref) it may either under or over-estimate the variance.} For time series, under appropriate conditions, a HAC estimator is consistent for the long-run variance of $\hat \psi_T$ and the averaged $\hat \psi^S_T(\theta_n^S)$ respectively. Computing the long-run variance for the averaged $\sum_{s=1}^S \psi(y_t^s,\dots,y_{t-L}^s)/S$ ensures that the estimate is consistent for both standard and antithetic draws. As before, an estimate for the non-averaged moment may not be consistent for antithetic draws because of the dependence between simulated moments.
For the Scrambled Method of Moments, the variance should not be computed as above because scrambled draws are not independent from one another. Theorem (ref) implies that the asymptotic variance only involves $\hat \psi_n$ in most cases; because simulation noise is asymptotically negligible.\footnote{See e.g. Proposition (ref).} One approach is to only compute the variance of $\hat \psi_n$. However, as illustrated in Section (ref), even though the simulation noise can be small in finite samples, it may not be completely negligible for some DGPs. In these cases, one would want to account for the variance attributable to $\hat \psi_n^S$. As discussed in Section (ref), to consistently estimates the variance of $\hat \psi_n^{S}$ one can evaluate $\hat \psi_n^{S}$ several times with different seeds for the scramble and compute the variance across these estimates.
In the following $\hat \theta_n^S$ and $\hat \theta_T^S$ will denote the scrambled estimator for static and dynamic models respectively. Consistency and asymptotic normality results are provided for the algorithms described above. The first set of assumptions below is standard in the Monte-Carlo simulation-based estimation literature.
To simplify notation, let: \[ \tilde \psi(x_i,\theta,u_i) \overset{def}{=} \psi( g(x_i,\theta,u_i),x_i ) \]
Assumption (ref) provides sufficient conditions to prove a uniform law of large numbers (ULLN) for $\tilde \psi$ and $\partial_\theta \tilde \psi$ using the scramble. The proof is similar to Jennrich1969.
Given the ULLN for the simulated moments and the fast convergence rate in Theorem (ref), the estimator is consistent and asymptotically normal. The main difference with standard SMM is that here the simulations do not inflate the asymptotic variance, even for $S=1$, whereas the simulation noise in SMM implies an additional $1/S$ factor.
As discussed in Section (ref), moments with covariates do not quite fit the setting described in Section (ref). Indeed, the scrambled draws are identically distributed but not independent. With the introduction of covariates, $\tilde \psi(x_i,u_i;\theta),i = 1,\dots,n$ are neither identically distributed nor independent which makes deriving ULLNs and CLTs challenging. Furthermore, if the moments are non-smooth in $\theta$ then the approach of Jennrich1969 cannot be applied and empirical process methods are required.
The main idea is to split the sample moments and the empirical process into two parts: one is non-identically distributed but independent and the other is identically distributed but not independent. The former can be handled using CTLs and empirical process results for heterogeneous arrays and assuming the later is smooth in $\theta$, it can be handled using the steps in Jennrich1969 as in Proposition (ref). The main assumption there is that integrating over $x_i$, while conditioning on $u_i$, transforms non-smooth into smooth moments. This puts restrictions on the moments and covariates used in the estimation.
Assumption (ref) i.-ii. ensure the Lindeberg condition holds for the heterogeneous array which is required to apply a CLT and the Jain-Markus Theorem VanderVaart1996. Conditions iii-iv. ensures that Theorem (ref) can be applied to the smoothed moments, i.e. after integrating out the covariates.
Proposition (ref) is similar to Pakes1989 with scrambled instead of MC draws. The variance $\tilde V$ can be computed using the steps described in Section (ref).
The following extends the results from Proposition (ref) to the Indirect Inference estimator of Gourieroux1993. The moments $\hat \psi_n,\hat \psi_n^S(\theta)$ are now defined as sample and simulated M-estimators:
Again, to simplify notation consider: \[ \tilde m(u_i,\theta;\psi) \overset{def}{=} m(y_i^s(\theta);\psi). \] As in Gourieroux1993, the binding function $\psi_{\infty}(.\cdot)$ is defined as: \[ \psi_{\infty}(\theta) \overset{def}{=} \text{argmin}_{\psi \in \Psi} \mathbb{E}\left[ \tilde m(u_i,\theta;\psi) \right]. \] Rather than matching sample moments, the indirect inference estimator matches functions (minimizers) of sample moment functions. Assumption (ref) below is more detailed than the high-level conditions in Gourieroux1993. Using implicit function arguments, it allows to express the estimator $\hat \theta_n^S$ in terms of the sample moments $\partial_\theta M_n^S$ which fit the setting of Section (ref) so that, eventually, Theorem (ref) applies.
Proposition (ref) is similar to the results found in Gourieroux1993 but here the simulation noise has no effect on the asymptotic variance as in Proposition (ref).
For simplicity, write: \[ \hat \psi_T^S(\theta) = \frac{1}{T S} \sum_{t=1}^{T S} \tilde \psi(u_t;\theta), \] where $u_t$ has the appropriate dimension $d$ given in Section (ref). For the qMC-only estimator, $\hat \psi_T^S$ is simply a cross-sectional average over short-time series. This fits the framework of Section (ref) directly and under the conditions in Assumption (ref) the estimator is consistent and asymptotically normal as shown in the Proposition below. As for the static case, the asymptotic variance is not inflated by the simulation noise, even for $S=1$.
For simplicity, write: \[ \hat \psi_T^S(\theta) = \frac{1}{T S} \sum_{t=1}^{T S} \tilde \psi(y^1_t,z^1_t,u_t;\theta), \] where $y_t^1,z_t^1$ are simulated using MC methods as in Algorithm (ref). The hybrid MC-qMC approach relies on MC simulations to approximately draw initial values from the ergodic distribution and is combined with the scramble to simulate a cross-section of paths.
Assumption (ref) i. is the usual geometric ergodicity condition duffie-singleton. Conditions ii.-v. are more restrictive, they hold if the moments are bounded. To relax these conditions, one would need to extend the CLT in Theorem 5.20 of White1984 to unbounded non-identically distributed dependent arrays which is outside the scope of this paper. Condition vi. requires the variance to be non-degenerate to apply a CLT. Otherwise, simulation noise is negligible in some directions which is not problematic in this setting.
Proposition (ref) is similar to duffie-singleton, the main idea to prove the result is to write the simulated moments as the sum of a mixing non-identically distributed heterogeneous array and an average of identically distributed non-independent terms. As in Proposition (ref), the former is handled using a specific CLT while the former uses Theorem (ref) with similar steps to Jennrich1969.
The following illustrates the finite sample properties of the Scrambled Method of Moments and Scrambled Indirect Inference in several simple examples and one application drawn from the heterogeneous agents literature. All simulations were carried out in R and C++ using the Rcpp package. Scrambled sequences were generated using the fOptions package.
The first example, drawn from Gourieroux1993, considers the estimation of a sample mean and variance of for an iid Gaussian sample: \[ y_i = \mu + \sigma e_i, \quad e_i \sim \mathcal{N}(0,1). \] This example illustrates Algorithms (ref), (ref) and Proposition (ref). As in the original paper, the auxiliary parameters $\hat \psi_n$ are the sample mean and variance of $(y_1,\dots,y_n)$: \[ \hat \psi_n = ( \hat \mu_n, \hat \sigma_n^2 )^\prime = \frac{1}{n} \sum_{i=1}^n (y_i, [y_i - \hat \mu_n]^2)^\prime. \] In the $5,000$ Monte-Carlo replications, the sample size is $n=100$ and $\theta_0=(\mu_0,\sigma_0^2) = (0,1)$. The number of simulated samples is $S=1,2,4$ and $20$. For SMM, $e_i^s$ is drawn using the random number generator rnorm in $R$ and antithetic draws are generated for $S=2,4$ and $20$ by taking $e_i^{s+S/2} = -e_i^{s}$ for each $s=1,\dots,S/2$. The fOptions package generates the scrambled Gaussian shocks directly. Table (ref) summarizes the biases and standard deviations of the estimators.
Because it has no simulation noise, the Method of Moments (MM) estimator has the smallest variance. SMM has a bias correction property for $\hat \sigma^2$ Gourieroux1993. For $\hat \mu_n$, antithetic draws and the scramble perform equally well for $S=2$ and $S=1$, respectively. For $\hat \sigma_n^2$, antithetic draws perform worse than SMM and the scramble. This is in line with the discussion in Section (ref). The scramble performs similarly to the MM while SMM requires $S=20$ to perform similarly. SMM and antithetic draws reduce the bias while the scramble does not. This reflects the fact that the scrambled $\hat \psi_n^S(\theta_0)$ approximates the asymptotic binding function $\psi_\infty(\theta_0) = \lim_{n \to \infty} \mathbb{E}(\hat \psi_n)$ while SMM and antithetic draws approximate the binding function $\psi(\theta_0) = \mathbb{E}(\hat \psi_n)$ which provides some finite sample bias correction.
The second example illustrates Algorithm (ref) with non-smooth moments and covariates. The DGP is a Probit model: \[ y_i = \mathbbm{1} \left\{ \theta_1 + \theta_2 x_i \theta + e_i \geq 0 \right\}, \quad e_i \overset{iid}{\sim} \mathcal{N}(0,1), \quad x_i \sim \mathcal{N}(0,1). \] The moments $\hat \psi_n$ consist of the intercept and the slope in an OLS regression of $y_i$ on $x_i$. In the $5,000$ Monte-Carlo replications, the sample size is $n=1,000$ and $\theta_0=(\theta_{1,0},\theta_{2,0}) = (1,1)$. The number of simulated samples is $S=1,2,4$ and $20$. The standard deviations of the estimators are reported in Table (ref).
The Scrambled Method of Moments outperforms SMM for $S=1$ and above. For $S \geq 2$, the scramble performs similarly to antithetic draws for estimating $\theta_0$ and $\theta_1$. The gains are less substantial than in the previous example.
To illustrate Algorithms (ref) and (ref) consider the following ARMA(1,1) model:\footnote{In the notation of Section (ref), the model can be written as: $y_t = \rho y_{t-1} + \sigma [z_{t,1} + \vartheta z_{t,2}],$ $(z_{t,1},z_{t,2})^\prime = (e_{t},z_{t-1,1})^\prime,$ with $z_{t} = (z_{t,1},z_{t,2})^\prime$.} \[ y_t = \rho y_{t-1} + \sigma[ e_t + \vartheta e_{t-1}], \quad e_t \overset{iid}{\sim} \mathcal{N}(0,1), \] In the $5,000$ Monte-Carlo replications, the sample size is $T=200$ and $\theta_0=(\vartheta_0,\rho_0,\sigma_0^2) = (0.5,0.5,1)$. The number of simulated samples is $S=1,2$. The moments are the OLS coefficients from regressing $y_t$ on its first $L=4$ lags and the variance of the OLS residuals. Using auto-covariances as moments instead yields similar results.
Algorithm (ref) requires sampling $(y^1_t,e^1_t)$ from its stationary distribution directly. The marginals are known since $e^1_t \sim \mathcal{N}(0,1)$ by assumption and $y^1_t \sim \mathcal{N}(0, [1+\vartheta^2 + 2\rho \vartheta]/[1-\rho^2]\sigma^2)$. Since they are jointly Gaussian, it is sufficient to compute their covariance, $\text{cov}(e_t^1,y_t^1) = \rho \vartheta \sigma$, to find their joint distribution: \[ \left(
\right) \sim \mathcal{N} \left( \left(
\right), \left(
\right)\right). \] Transforming independent bivariate scrambled Gaussian shocks into draws from the joint distribution above is then straightforward. For Algorithm (ref), the $(y^1_t,e^1_t)$ need to be sampled using MC methods. First, the initial value $(y^1_0,e^1_0)=(0,0)$ is set and a path $(y_t^1,e_t^1)$ is simulated with random MC draws. Once these $(y_t^1,e_t^1)$ are simulated, the remaining $(y_t^2,e_t^2,\dots,y_t^5,e_t^5)$ are computed using scrambled Gaussian shocks.\footnote{Here the number of lags used to compute the moments is $L=5$ because $y_t^5$ is regressed on its 4 lags $y_t^4,\dots,y_t^1$.}
Table (ref) compares the Maximum Likelihood Estimator (MLE) with SMM and antithetic draws, the qMC-only scramble from Algorithm (ref) (reported in the Scramble column) and the hybrid MC-qMC scrambled from Algorithm (ref) (reported in the Scramble-MC column).
MLE corresponds to the lower bound for variance of the estimators. The qMC-only scramble from Algorithm (ref) outperforms SMM and antithetic draws. Antithetic draws perform worse than SMM using the same $S$ which further illustrates the discussion in Section (ref). The hybrid MC-qMC Algorithm (ref) performs better than SMM and antithetic draws and, as expected, worse than the qMC-only approach.
The last example is a more substantial model borrowed from BROWNING2010.\footnote{The data generating process considered here involves all the coefficients found in BROWNING2010, Table 2 minus the measurement errors and the time-trend in the ARCH component which are not considered in this Monte-Carlo exercise.} Simulation-based estimation is commonly used in this heterogeneous agents literature due to the complexity and intractability of the models.\footnote{See e.g. Guvenen2011 for an overview of the computation and estimation of heterogeneous agents models.} The baseline data generating process is an ARMA(1,1) at the individual level:
where the drift $\alpha_i$, long-run mean $\delta_i$, AR and MA coefficients $\beta_i,\theta_i$ as well as the persistence coefficient $\omega_i$ all vary at the individual level. The Gaussian shocks to log-income in the time dimension are denoted by $\varepsilon_i$ while the Gaussian shocks to the ARMA coefficients $\alpha_i,\beta_i,\dots$ will be denoted by $\eta_i$. The initial value for log-income $y_{i,0}$ is drawn as:
The heterogenous ARMA coefficients are then drawn using:
where $\text{logit}$ is the usual logistic transformation $\text{logit}(x) = 1/[1+\exp(-x)]$. $\eta_{i,0},\dots,\eta_{i,2} \overset{iid}{\sim} \mathcal{N}(0,1)$. For a discussion of the parameters and the role of the transformations, see BROWNING2010. $\nu_{i,0}$ is the initial value for the ARCH-type heteroskedasticity in the shocks $\varepsilon_{i,t}$:
where $e_{i,0},\dots,e_{i,T} \overset{iid}{\sim} \mathcal{N}(0,1)$. In the simulations, the number of households is $n=1,000$; the number of time periods is $T=30$. As in the original paper, a burn-in period of $T_{burn}=3$ periods is used to reduce the effect of the initial conditions. The parameter values are taken from Table 2 in BROWNING2010 and the moments are those described in their Appendix A.2 except the ones involving year of birth which are not considered in these simulations. In a nutshell, the moments involve the aggregation of individual-level OLS coefficients, moments based on OLS residuals, autocorrelations and measures of social mobility.
The implementation of SMM is standard and described in Appendix A.4 of the original paper. For the scramble, a $(n\times S) \times (T+T_{burn} + 3) = (1,000 \times S) \times 36$ matrix of scrambled standard gaussian shocks is drawn. The integration dimension $d=36$ is sufficiently large to illustrate the finite sample performance of the scrambled method of moments with a relatively large number of shocks. The first three dimensions (columns of the matrix) correspond to $\eta_{i,0},\dots,\eta_{i,2}$, the remaining dimensions correspond to time dimensions $e_{i,1},\dots,e_{i,T+T_{burn}}$. The rows correspond to the cross-sectional dimension of the shocks, i.e. the $i=1,\dots,n\times S$ index.
The results from the $2,000$ Monte-Carlo replications are presented in Table (ref) for $S=1,2,4$. SMM and antithetic draws are used as a benchmark for the scramble with either a large sample of $n \times S$ individuals (as in Algorithm (ref) or $S$ samples of $n$ individuals (as in Algorithm (ref)). The scramble generally outperform SMM and antithetic draws. Antithetic draws either under or over-performs SMM depending on the parameter of interest which is in line with previous discussions. Both implementations of the Scrambled Method of Moments perform similarly. For some coefficients, there is little to no improvement in increasing $S$ from $1$ to $2$ or $4$. For most coefficients, the scramble with $S=2$ outperforms SMM with $S=4$. Furthermore, using the same $S=4$, the replications were computed about $15\%$ faster for the scramble than SMM. Since the only difference between the two is the shocks used in the simulations, this reflects faster convergence of the optimizer. Possibly because the scramble are smoother (less noisy) than the MC moments which makes the objective function easier to minimize.
This paper proposes several algorithms implementing Owen's scramble for simulation-based estimation. Since the method is designed for computing integrals of iid sequences, some care is needed when simulating data with covariates or time series. Large sample results are provided to support the proposed algorithms. The results for dynamic models could be extended to non-smooth bounded moments through additional stochastic equicontinuity results using the inequality in Andrews1994 for instance. The simulations illustrate the finite performance of the Scrambled Methods of Moments and Scrambled Indirect Inference compared to other commonly used methods. The last example suggests the scramble could be useful in larger scale problems found in the heterogenous agents literature where SMM is commonly used.