Problems using randn in Matlab 2010 to simulate an SDE path path

[Lorenzo]

Hello,

I must premit I am not a numerical analyst/Matlab expert so please forgive me if I do not express myself in the appropriate terms.

Few months ago I started working on a semi-analytical approximations to an expectation of a functional of a terminal SDE value (yes, financial maths). I tested them against a Monte Carlo benchmark built upon a Euler scheme for a particular SDE and arithmetic means of the path-ways terminal values. I used Matlab 2007, do not remember the version though. I came up with 3 different analytical formulas approximating the Monte Carlo simulation (100000 path samples) exactly to the second-third decimal point.

I have now changed to the 2010 Matlab version. Running the same .m files for the Monte Carlo simulati now yields a 0.04-0.1 error with the semi analytical approximation. In relative terms the new expectations have an average 1% error respect to the old ones, systematically. This makes me suspect the algorithm I am using combined with the 2010 randn algorithm introduces a bias.

I have tried the following:

1) Reseeding the stream, as explained in the math help, via

RandStream.setDefaultStream(s);
s = RandStream('mt19937ar','Seed', 5489);

after starting the session, at global level before the path simulations, after each path simulation, after each gaussian sample for the sde. The first had no effect, the latter three killed the randomness.

2) Defining the algorithm in both of the following ways.

2.1) Generating a global random matrix on which the SDE process operates.
2.2) Generating one normal sample for each mesh of the path of the process. No difference.

As a remark my original code implemented on my old 2007 version used the second ( I know it is not optimal) with no initial seeding for randn. I didn't even consider polishing it becuase results were already strikingly accurate and runtime was not an issue.

I would appreciate any help or suggestion. I know it is a bad thing to say I do not really know how random number generators work, I use them as a black box. I guess I am learning something from this situation.

Lorenzo

[Peter Perkins]

I have now changed to the 2010 Matlab version. Running the same .m files
for the Monte Carlo simulati now yields a 0.04-0.1 error with the semi
analytical approximation. In relative terms the new expectations have an
average 1% error respect to the old ones, systematically. This makes me
suspect the algorithm I am using combined with the 2010 randn algorithm
introduces a bias.

It's possible, but there are other simpler explanations too.

1) Reseeding the stream, as explained in the math help, via
RandStream.setDefaultStream(s);
s = RandStream('mt19937ar','Seed', 5489);
after starting the session, at global level before the path simulations,
after each path simulation, after each gaussian sample for the sde. The
first had no effect, the latter three killed the randomness.

The order of those two lines should be reversed, but modulo that, you
would expect to see what you did see, right?

a) Are your results stable with respect to multiple runs of the
simulation? What happens if you run your simulation several times
without touching anything to do with the random number generator?


b) What happens if you try another one of the new (since R2008b) generators:

s = RandStream('mrg32k3a','Seed', 0);
RandStream.setDefaultStream(s);


c) Or if you use the RNG that was underneath RANDN in R2007a:

s = RandStream('shr3cong','Seed', 0);
RandStream.setDefaultStream(s);


c) And of course, if you do something like

randn('state',0)

you will get exactly the old (pre-R2008b) behavior of MATLAB's random
number generation, so if this is a RNG issue you should see your old
results.

Hope this helps.

[Lorenzo]

Peter Perkins <***@MathRemoveThisWorks.com> wrote in message <ice0b6$rod$***@fred.mathworks.com>...

Thanks for your reply Peter, I much appreciate it

Post by Peter Perkins
a) Are your results stable with respect to multiple runs of the
simulation? What happens if you run your simulation several times
without touching anything to do with the random number generator?

Yes, and that precisely what worries me. To give you an example on Matlab 7.3 a result was stable to

6.714 + err~10^-4

On Matlab 2010 and EVEN 7.5 (therefore pre 2008 RNG) the same simulation, ran multiple times, with the same number of samples, is stable to

6.677+ err~10^-4

Post by Peter Perkins
s = RandStream('mrg32k3a','Seed', 0);
RandStream.setDefaultStream(s);
s = RandStream('shr3cong','Seed', 0);
RandStream.setDefaultStream(s);
c) And of course, if you do something like
randn('state',0)

I tried all of this and unfortunately nothing would work.

One thing I forgot to mention. My older simulations were performed on a 2007 STUDENT VERSION. To your (or anyone else that might be reading) knowledge, are there any reported unfixed bugs/problems for such a version, in particular related to RNGs?
Even so, assuming that there are, I have heaps of papers and .m files of maths showing that, for instance in the example above, the correct value is

6.714 +Numerical approx error+Analytical approx error~10^-4

and not

6.667+ Numerical approx error+Analytical approx error~10^-4

I am obviously feeling very uneasy.
I might be going back ato the older version and try to fiddle around with randomness THERE to see if there is any sensitivity to change of seeding/randomness algorithms.

Lorenzo

[Peter Perkins]

Post by Lorenzo

Post by Peter Perkins
s = RandStream('shr3cong','Seed', 0);
RandStream.setDefaultStream(s);
c) And of course, if you do something like
randn('state',0)

I tried all of this and unfortunately nothing would work.

The default behavior for random numbers in MATLAB changed twice since
MATLAB Version 5. In R2007a, the default uniform generator underneath
RAND became the Mersenne Twister, but the normal generator underneath
RANDN remained unchanged. In R2008b, RANDN stopped using its own
separate generator and now uses the ziggurat transformation algorithm on
top of the same Mersenne Twister as RAND.

I'm not 100% sure where your Student Version falls, but if you see this
in a new session

Post by Lorenzo

Post by Peter Perkins
rand

ans =
0.81472

then it's the MT, if you see 0.95013, it's the earlier one. If you
execute this

rand('twister',0) % to mimic R2007a-R2008a
randn('state',0)

or this

rand('state',0) % to mimic R2006b and earlier
randn('state',0)

the values coming from RAND and RANDN will be identical to those in
earlier versions. If you do that, and still see differences from when
you ran under your Student Version, then the problem is not in the
random number generator. You should be able to do the above and use the
debugger to step through one iteration of your MC simulation in both
versions, and find out exactly where the differences are.

Hope this helps.

[Lorenzo]

Peter Perkins <***@MathRemoveThisWorks.com> wrote in message <ice0b6$rod$***@fred.mathworks.com>...

Thanks for your reply Peter, I much appreciate it

Post by Peter Perkins
a) Are your results stable with respect to multiple runs of the
simulation? What happens if you run your simulation several times
without touching anything to do with the random number generator?

Yes, and that precisely what worries me. To give you an example on Matlab 7.3 a result was stable to

6.714 + err~10^-4

On Matlab 2010 and EVEN 7.5 (therefore pre 2008 RNG) the same simulation, ran multiple times, with the same number of samples, is stable to

6.677+ err~10^-4

Post by Peter Perkins
s = RandStream('mrg32k3a','Seed', 0);
RandStream.setDefaultStream(s);
s = RandStream('shr3cong','Seed', 0);
RandStream.setDefaultStream(s);
c) And of course, if you do something like
randn('state',0)

I tried all of this and unfortunately nothing would work.

One thing I forgot to mention. My older simulations were performed on a 2007 STUDENT VERSION. To your (or anyone else that might be reading) knowledge, are there any reported unfixed bugs/problems for such a version, in particular related to RNGs?
Even so, assuming that there are, I have heaps of papers and .m files of maths showing that, for instance in the example above, the correct value is

6.714 +Numerical approx error+Analytical approx error~10^-4

and not

6.667+ Numerical approx error+Analytical approx error~10^-4

I am obviously feeling very uneasy.
I might be going back ato the older version and try to fiddle around with randomness THERE to see if there is any sensitivity to change of seeding/randomness algorithms.

Lorenzo