Coding of Ito Stochastic Process

Question

I am trying to implement a routine in mathematica/matlab for a stochastic process. Any code written here is for mathematica, but if someone can help me with encoding this in matlab (if they're more familiar with that) then that would be fine as well. However, mathematica is the priority if possible.

I will state what I want to get after going over the equations first.

The following is the Ito stochastic Process that I am interested in (where z(t)=[x(t),y(t)]:

with the following quantities:

We also have the following (from now I will type x(t),y(t) as x,y):

where, is the first exit time and GT(x) is a smooth function of x (please see at the end)

Goal: I want to get Q0 in terms of x and only.

=============================================================================== To find the first exit time maybe the following can be used (courtesy of b.gatessucks). Please note that the following is a mathematica code.

CONSTRAINT FOR EXIT TIME:

const[x_, y_] := And[10^-8 <= y <= 10^-3, 0.9*(Uc) <= a1*x^2/2 - a3*x^4/4 <= 1.1*(Uc)]

x0 = 0.1;   (* starting point for x[t] *) y0 = 0.1; (* starting point for y[t] *)

proc = ItoProcess[  {\[DifferentialD]x[t] == y[t] \[DifferentialD]t,   \[DifferentialD]y[t] == (-G*y[t] - (a1*x[t] - a3*x[t]^3) - eps*b3b*y[t]^3) \[DifferentialD]t + Sqrt[2*eps*G] \[DifferentialD]w[t]},    {t, x[t], y[t], Boole[const[x[t], y[t]]]},   {{x, y}, {x0, y0}}, {t, 0}, w \[Distributed] WienerProcess[]]

Exit time is found:

SeedRandom[3] 
sim = RandomFunction[proc, {0, 1, 0.001}];
First@Select[sim[[2, 1, 1]], #[[4]] == 0 &]

that outputs {t, x[t], y[t], Boole[const[x[t], y[t]]]}

I need to be able to use the above code to find the exit time and then use it in Q0. The above snippet to find exit time produces non zero times for properly chosen .

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

Numerical Task that needs to be set up to find Q0(x):

--> We start from the integral term in this expression. First find 's for many initial conditions (say 100 for now -- i.e., 100 exit times) starting from inside the domain (maybe by using the snippet already mentioned in this post). Now the integral can be evaluated for each of the 100 exit times as functions of x and .

--> Now, the expectation of the integral in Q0 is a sample average (by Law of Large Numbers) of all the 100 integrals evaluated before. Thus, Q0 can be found and it should only be a function of x and .

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

My issues:

1. The code above for the exit time only seems to produce non-zero exit times for appropriately chosen initial conditions. If someone can elucidate as to how to pick the appropriate such that there will be sufficient exit time being produced then that would be appreciative. I really want to keep my constraints as specified above in const[x_, y_], but if there seems to be no hope in finding tractable results then I wouldn't mind relaxing it.

2. The code below for GT(x) results in singularities -- I received error message from DSolve regarding indeterminate expressions....both b0[xc] and DrhoDy[xc] become indeterminate and so the DSolve gives problems when using the initial condition GT[xc] as it becomes indeterminate as well...any way around this matter would be gladly appreciated.

3. Finally, I really need someone's help in evaluating the 100 integrals in mathematica efficiently (since the terms are huge) for each of the exit times found previously and taking the expectation. I am not sure as to how to find the exit times properly.

================================================================================ To find GT(x):

b1b = 0.9; b3b = .8; a1b = 0.1; a3b = 0.2; G = (1/0.1^2)*
  b1b; a1 = (1/.1^2)*a1b; a3 = (1/.1^2)*a3b; xc = Sqrt[a1/a3]; Uc = 
 a1*xc^2/2 - a3*xc^4/4; L = (G + Sqrt[G^2 + 4*(-(a1 - 3*a3*xc^2))])/2;
y[x_] := (x *a1 - x^3*a3)/(G + L)

b0[x_] := (y[
      x]*(a1*x \[Minus] 
       a3*x ^3)*(1 \[Minus] (a1 \[Minus] 3*a3*x ^2)) \[Minus] 
    G*y[x]*(a1 \[Minus] 3*a3*x ^2)) /(y[
      x] ^2 + (\[Minus]G*y[x] + a1*x \[Minus] a3*x ^3) ^2)


DrhoDy [x_] := 
 Sqrt[y[x]^ 2 /(y[x] ^2 + (\[Minus]G*y[x] + a1*x \[Minus] a3*x ^3) ^2)]


linearequation = 
  y[x]*GT '[x] + b0[x]*GT[x] == G*DrhoDy [x]^2 *GT[x]^ 3;

GT[x] = 
 DSolve[{linearequation, GT[xc] == Sqrt[b0[xc]/( G*DrhoDy [xc]^2 )]}, 
  GT, x]

ODE:

that satisfies the Initial Condition:

where,

and