# Quasi Monte Carlo in Matlab

I want to use Quasi Monte Carlo to try and improve the convergence of a simulation I am running.

The random numbers are simply to produce the observation errors for a standard linear regression model. Which is then estimated using a number of different regression techniques. This is done repeatedly to estimate the mean square error of each model.

I'm fairly new to Quasi Monte Carlo but is is likely to help in this situation I am just using it to produce 10k random numbers. It seems that generally I can expect quicker convergence of the order of (1/n) rather than n^(-0.5):

http://en.wikipedia.org/wiki/Quasi-Monte_Carlo_method

However it also states that the QMC numbers are not truly random, so I just wonder what the implications might be for any statistical tests I might want to run on the results.

1.) I guess what I want to know are the pros and cons of MC v QMC. (would you always want to use QMC if its available?) 2.) What tests can I use to ascertain which is best for my application? (seems any test that depends on the numbers being truly random will fail?)

I know that this can be done in Matlab using

q = qrandstream('halton',NSteps,'Skip',1e3,'Leap',1e2); RandMat = qrand(q,NRepl); z_RandMat = norminv(RandMat,0,1);

which is taken from this paper.

http://papers.ssrn.com/sol3/cf_dev/AbsByAuth.cfm?per_id=1519239

It seems there are other low discrepancy numbers such as Sobol sequence available in Matlab and again would just like to know what tests I can use to ascertain which is best for my situation.

• Could you explain a bit more what you are exactly trying to achieve? E.g. what do you mean by "convergence"? It is not clear to me what is supposed to converge to what. The second paragraph sounds like you were using the Quasi random numbers as "noise" which you add to something so that your regression methods have something to pick up. The convergence results you quote are related to the variance of the standard MC-Estimator which is something different from regression. – g g Jan 27 '14 at 12:54
• I am using QMC to generate different realizations of the noise parameter to be used in a regression model. By reestimating the model repeatedly based on these different realisations of the noise I can build up an estimate of the mean square error of the estimator. Your point about the difference between standard MC-estimator convergence and what I am trying to is really the crux of the matter. Also I'm trying to this on just 5000 iterations if poss, which I know is kinda borderline anyway. From the results I am seeing it seems to be slightly better but would like to be able to quantify this. – Bazman Jan 27 '14 at 16:28

One of the main things you give up is a simple halting condition for your estimation algorithm. With pseudorandom numbers, the algorithm can keep track of the standard error, and stop when it has passed a threshold:

error_est = Inf
n = 0
while not error_est < target_precision:
n = n + 1
x = new_random_sample()
samples.append( F(x) )
error_est = 3 * std_deviation( samples ) / sqrt(n)
value_est = mean( samples )


Since with QR sequences your error estimate can't be taken from the standard deviation, you cannot apply this algorithm. Instead, practitioners often either just choose a number of samples to take a priori, or set up a halting condition that checks for the running mean to be changing very little.

As for tests of one sequence versus another, the typical approach is to set up a large suite of sample problems. The winner is then the QR sequence that has the best computation time to precision tradeoff on the problem suite.

When we did this for equity exotics in the 1990s, we found Niederreiter sequences beat Sobol and Faure, though by a very small amount.

• Sorry to be thick but why can't you estimate the error (I assume you mean standard error) from the standard deviation? I was hoping to compare the different schemes: pseudorandom, Hald and Sobol by looking at the standard deviations of the estimators. Does your comment above mean this approach is invalid? It appears the primary reason to use QMC is to reduce the standard deviation of the estimates so there must be some way of doing it? – Bazman Jan 24 '14 at 15:36
• The standard deviation of the estimator is a different beast from the standard deviation of the samples. You can indeed get an idea of the estimator's error by running on subsamples. However, unlike with PR sequences, you will find the estimator errors do not scale as $N^{-1/2}$. – Brian B Jan 24 '14 at 17:31
• The QMC variables scale 1/N (thoeretically) on subsamples. Are you saying that this doesn't necessarily mean that the standard of the deviation will scale in the same way? If so how does it scale? – Bazman Jan 24 '14 at 19:22
• A sequence of pseudo-random numbers are i.i.d. samples, and hence you can apply the central limit theorem to estimate the variance (error). However, for a quasi-random sequence the samples are not independent, and hence one cannot apply the CLT. But if you take one sequence and randomise/scramble it, the estimates from these scrambled sequences do form i.i.d estimates, and can be combined to form an error estimate. – oliversm Aug 23 '18 at 22:30
• @oliversm is not correct, except in the limit as $N \rightarrow \infty$. The entire point of quasirandom sequences is that they cover the space differently. – Brian B Aug 24 '18 at 23:47