Trying to simulate a birth death process Question:
A birth and death process is a continuous time
Markov chain. 
Find an approximative numerical value for the probability 
P {max0≤t≤10 X(t) ≥ 10} 
for a birth and death process {X(t)}t≥0 
with birth intensities λ0 = λ1 = λ2 = . . . = 1
and death intensities µ1 = µ2 = µ3 = . . . = 2 
that is in steady-state (that is, which is
started according to its stationary distribution).
My attempt at a solution:
rep = 100000

count = 0


for (i in 1:rep) {

   time = 0
   xt = 0
   succ = 0
   start = 0

   u.rand = runif(1)

   # These if statements determine where we start 
   # according to the stationary distribution

   if (u.rand<1/2) {start = 0}                   

   if (1/2<=u.rand && u.rand<3/4) {start = 1}

   if (3/4<=u.rand && u.rand<7/8) {start = 2}

   if (7/8<=u.rand && u.rand<15/16) {start = 3}

   if (15/16<=u.rand && u.rand<31/32) {start = 4}

   if (31/32<=u.rand && u.rand<63/64) {start = 5}

   if (63/64<=u.rand && u.rand<127/128) {start = 6}

   if (127/128<=u.rand && u.rand<255/256) {start = 7}

   if (255/256<=u.rand && u.rand<511/512) {start = 8}

   if (511/512<=u.rand && u.rand<1023/1024) {start = 9}

   if (1023/1024<=u.rand) {start = 10}

   # xt is the position during the while loop
   # and its a success if xt gets bigger than 10 

   xt = start                                          

   while ((time<=10) && (succ==0)) {

    mu = rexp(1, 2)                                   # death rate
    la = rexp(1, 1)                                   # birth rate


    # Here we get a birth and xt increases by 1.
    if (la<mu) {                                      

    xt = xt + 1

    time = time + la

    # Here we get a death as long as xt not is equal to 0

    } else if ((mu<la) && (xt!=0)) {                 

    xt = xt - 1

    time = time + mu

    # at xt=0 we can only get births

     } else if (xt==0) {                               

     xt = xt + 1

     time = time + la }


      # Here we register the successes.

      if (xt>=10) {succ = 1}                              

}

count = count + succ                  # And here we add them up

}

count/rep                                             

# And this frequency should be equal to
# the sought after probability

I get count/rep = 0.00615 but it should be 0.0084865
Something is wrong in my code.
Can anyone help me out? 
Im pretty sure im doing something wrong in the while-loop.
 A: First of all, are you sure that 0.0085 is the correct answer? Looking at here suggests otherwise???
Some comments about your code:
(1) I think your code is essentially doing the right thing. However, the structure of the code within your while loop makes it difficult to read. I would split up the code into two 'if' clauses, i.e. if(xt==0) {...} and if(xt!=0) {...}, reflecting the special nature of the xt==0 situation, so that it is clear what actions are performed for these two different cases.
(2) Be careful as it is possible that 'time' may exceed 10 in your while loop, so you would need to ask that time<=10 as well as xt>=10 in order to register a success.
(3) u.rand>=1023/1024 corresponds to the situation where start>=10. Surely in these cases success is already guaranteed i.e. you should set succ=1.
(4) The statements about u.rand could be combined into a for loop, which would save typing time and reduce the risk of typing error.
(5) Running your code more than once reveals substantial between-run variability. You could reduce this by increasing 'rep'. Even better, you could stratify your sampling of the stationary distribution, as at the moment your use of 'runif' is introducing additional variability into your answers. To stratify, you could start by specifying rep=1024 (for example), then for the first 512 reps set start=0, for the next 256 reps set start=1, and so on.
