# Main simulation function

# simulates diffusion process of the form
# dX(t) = a(X(t),t)dt + b(X(t),t)dZ(t)
# where Z(t) is a standard Wiener process

Itosim_function (t = 50, dt = 0.01, afun = a, bfun = b, x0 = 1) 
{
            n <- t/dt
            xvec <- vector("numeric", n)
            x <- x0
            sdt <- sqrt(dt)
            for (i in 1:n) {
                    t <- i * dt
                    x <- x + afun(x, t) * dt + bfun(x, t) * rnorm(1, 0, 
                            sdt)
                    xvec[i] <- x
            }
            ts(xvec, deltat = dt)
}


# Examples of use
op_par(mfrow=c(4,2))

a_function(x,t){ 0.05 }
b_function(x,t){ 0.5 }
x_Itosim()
plot(x,main="Generalised Wiener process")
plot(exp(x),main="exponential of Generalised Wiener process")

a_function(x,t){ (0.05+0.5*0.5*0.5)*x }
b_function(x,t){ 0.5*x }
x_Itosim()
plot(x,main="Geometric generalised Wiener process")
plot(log(x),main="log of geometric generalised Wiener process")

a_function(x,t){ -3*(x-1) }
b_function(x,t){ 1 }
x_Itosim()
plot(x,main="Gaussian OU process")
plot(exp(x),main="exponential of Gaussian OU process")

a_function(x,t){ (0.5 + 3*1 - 3*log(x)) * x }
b_function(x,t){ x }
x_Itosim()
plot(x,main="Geometric Gaussian OU process")
plot(log(x),main="log of Geometric Gaussian OU process")

par(op)