Skip to content

Problems Simulating (PR#3471)

2 messages · smiles@dept.econ.yorku.ca, Peter Dalgaard

#
Full_Name: Stan Miles
Version: 1.5.0.
OS: Windows 2000 Server
Submission from: (NULL) (130.63.74.220)


Short Description of the Problem:  

I use R to simulate a random variable, with mean 0.14 and variance 0.2 .  I
simulate 20 sets, 30000 realizations/set.
I take the average of each set, and all 20 of the averages are 
much higher than 0.14.  In fact, they are All about 2-3 stdev higher. 
It seems that R is not generating the random variable with the mean I ask for.

The Details:

I simulate Prices [Pt] and Returns that are defined as 
Rt = [P(t)-P(t-1)]/P(t-1).  The returns have Normal distribution
=>  Rt is has N(mu=0.14, sigmaSquared=0.2) distribution.

To do the task above, I use the following files:

c:\ProgramFiles\R\rw1050\Code1.txt

consists of the following text:

====================
gpsim1_function(t,runs,mu,delta,stdev){
jold=0
xx_rep(0,t*runs)
for(i in 1:runs){
ind=(i-1)*t
xx[ind+1]=1
for(j in 2:t) {
xx[ind+j]_xx[ind+j-1]*
exp(rnorm(n=1,mean=mu*delta*delta,sd=stdev*delta))
}
}
xx

}

=======================

File "experiment4commands.txt" has text:

source("code1.txt")
aa=gpsim1(30000,1,0.14,0.052342392,0.4472135955)
sink("experiment4mu14z20Path1.txt")
for(i in 1:length(aa)){print(aa[i])}
sink()

source("code1.txt")
aa=gpsim1(30000,1,0.14,0.052342392,0.4472135955)
sink("experiment4mu14z20Path2.txt")
for(i in 1:length(aa)){print(aa[i])}
sink()

:
:
:

source("code1.txt")
aa=gpsim1(30000,1,0.14,0.052342392,0.4472135955)
sink("experiment4mu14z20Path20.txt")
for(i in 1:length(aa)){print(aa[i])}
sink()
=========================
{Above, 0.052342392 = 1/365, it is my discretization factor.}


I paste the text from the "experiment4commands.txt" file into R.

I now open each of the resulting 20 files in Excel.  The output [with 
variable Pt] is located in column B.  To get back Rt, I paste
the following formula into cells C2:C30000 - 
Ci = [B(i)-B(i-1)]/B(i-1).

My discretization factor was 1/365, thus, when I take the Average 
[and the Variance] of cells C2:C30000, I also have to multiply this 
number by 365, to get back the annual Rt.

I then made a table of the resulting 365*Average(C2:C30000), for each
of the 20 Sets that I have simulated using the methodology described above.


Set     365*Average(C2:C30000) for the set

1	0.284948
2	0.238827892
3	0.176492
4	0.242007
5	0.263882
6	0.211672
7	0.203536
8	0.19211
9	0.225185
10	0.234478
11	0.272792
12	0.205158
13	0.189349
14	0.222317
15	0.220652
16	0.165553
17	0.294211
18	0.333488
19	0.147971
20	0.290742

The average of the realizations of mu [true value = 0.14] over 
20 paths, each with 30000 observations, is about 0.23.

For this case, daily mu = mu over each algorithm step
is 0.14/365 = 0.0003835616438

The daily variance is 0.2/365 = 0.0005479452055

However, according to statistics, as the number of observations increases
variance will fall.  In particular, for 30,000 observations, the 
StDev is calculated as follows:  

SQRT(0.2/365)/SQRT(30000) = 0.0001351474 = {1 STDEV}

The average of the values in the table above, 0.23, is

[(0.23/365) - (0.14/365)]/0.0001351474 = 1.824 StDeviations above the mean.
As you can see Most of the 20 sets had realized values about 2-3 StDeviations
higher than the mean.

Something is wrong.

Thank you for considering this bug report.
#
smiles@dept.econ.yorku.ca writes:
This is not a bug report. It is a request for finding the bug or
logical error in your code. If you had had a real bug to report, you
should have investigated more carefully and boiled it down to a much
simpler example. In particular, it is not a good idea to require the
use of commercial software on one particular computing platform to
reproduce your results when the same operations could just as easily
have been done in R.

Your R version is getting a bit old, BTW. The pre-1.7.x random
number generators had problems, but not of this magnitude.