########################
###### EXERCISE 1 ######
########################
#### DETERMINISTIC #####
########################

## Step size:
dt   <- 0.1

## Define short time sequence:
time <- seq(0,5,dt)
N    <- length(time)

## Function returning the trajectory of X(t) as a function of
## the initial condition and the time points:
X.fun.det <- function(time.vector,X0){
  sign(X0)*exp(time.vector)*(X0^{-2}+exp(2*time.vector)-1)^{-0.5}
}

## Plot 4 trajectories against time for 4 different initial conditions:
plot(time,X.fun.det(time,X0 = 2),type="l",ylim=c(-2,2),ylab="X(t)")
lines(time,X.fun.det(time,X0 = 0.2),type="l")
lines(time,X.fun.det(time,X0 = -2),type="l")
lines(time,X.fun.det(time,X0 = -0.2),type="l")

########################
###### EXERCISE 2 ######
########################
# DOUBLE WELL POTENTIAL#
########################

## The potential for -2 < x < 2 :
x.seq <- seq(-2,2,0.1)
plot(x.seq,x.seq^4/4-x.seq^2/2,type="l",ylab="U(x)",xlab="x")

########################
###### EXERCISE 3 ######
########################
#### ADDITIVE NOISE ####
########################

## Function returning a trajectory of X(t) as a function of
## the initial condition and the size of the noise:
X.fun.stoch.add <- function(sigma,X0){
  X     <- 0*(1:N)
  X[1]  <- X0
  for(i in 2:N){
    dWt   <- sqrt(dt)*rnorm(1,0,1)    ## Draws a random increment
    X[i]  <- X[i-1] + (-X[i-1]^3 + X[i-1])*dt + sigma*dWt
  }
  X
}

## Plot 4 trajectories against time for 4 different initial conditions, each
## in a different panel, for the given value of the size of the noise, nysigma:
nysigma <- 0.5
par(mfrow=c(2,2))     ## Allows multiple panels in the same plot
plot(time,X.fun.stoch.add(sigma=nysigma,X0 = 2),type="l",ylim=c(-2,2),ylab="X(t)")
plot(time,X.fun.stoch.add(sigma=nysigma,X0 = 1/2),type="l",ylim=c(-2,2),ylab="X(t)")
plot(time,X.fun.stoch.add(sigma=nysigma,X0 = -1/2),type="l",ylim=c(-2,2),ylab="X(t)")
plot(time,X.fun.stoch.add(sigma=nysigma,X0 = -2),type="l",ylim=c(-2,2),ylab="X(t)")

## Same thing:
nysigma <- 0.1
par(mfrow=c(2,2))     ## Allows multiple panels in the same plot
plot(time,X.fun.stoch.add(sigma=nysigma,X0 = 2),type="l",ylim=c(-2,2),ylab="X(t)")
plot(time,X.fun.stoch.add(sigma=nysigma,X0 = 1/2),type="l",ylim=c(-2,2),ylab="X(t)")
plot(time,X.fun.stoch.add(sigma=nysigma,X0 = -1/2),type="l",ylim=c(-2,2),ylab="X(t)")
plot(time,X.fun.stoch.add(sigma=nysigma,X0 = -2),type="l",ylim=c(-2,2),ylab="X(t)")

## Same thing:
nysigma <- 1
par(mfrow=c(2,2))     ## Allows multiple panels in the same plot
plot(time,X.fun.stoch.add(sigma=nysigma,X0 = 2),type="l",ylim=c(-2,2),ylab="X(t)")
plot(time,X.fun.stoch.add(sigma=nysigma,X0 = 1/2),type="l",ylim=c(-2,2),ylab="X(t)")
plot(time,X.fun.stoch.add(sigma=nysigma,X0 = -1/2),type="l",ylim=c(-2,2),ylab="X(t)")
plot(time,X.fun.stoch.add(sigma=nysigma,X0 = -2),type="l",ylim=c(-2,2),ylab="X(t)")

########################
###### EXERCISE 4 ######
########################
# MULTIPLICATIVE NOISE #
########################

## Euler scheme ##
## Function returning a trajectory of X(t) as a function of
## the initial condition and the size of the noise:
X.fun.stoch.mul.eu <- function(sigma,X0){
    X     <- 0*(1:N)
    X[1]  <- X0
    for(i in 2:N){
      dWt  <- sqrt(dt)*rnorm(1,0,1)    ## Draws a random increment
      X[i] <- X[i-1] + (-X[i-1]^3 + X[i-1])*dt + sigma*sqrt(1+X[i-1]^2)*dWt
    }
    X
}

## Milstein scheme ##
## Function returning a trajectory of X(t) as a function of
## the initial condition and the size of the noise:
X.fun.stoch.mul.mil <- function(sigma,X0){
    X     <- 0*(1:N)
    X[1]  <- X0
    for(i in 2:N){
      dWt  <- sqrt(dt)*rnorm(1,0,1)    ## Draws a random increment
      X[i] <- X[i-1] + (-X[i-1]^3 + X[i-1])*dt + sigma*sqrt(1+X[i-1]^2)*dWt +
              0.5*sigma^2*X[i-1]*(dWt^2-dt)
    }
    X
}

nysigma <- 0.5

## Plot 4 trajectories against time for 4 different initial conditions, each
## in a different panel, simulated with the Euler scheme:
par(mfrow=c(2,2))     ## Allows multiple panels in the same plot
plot(time,X.fun.stoch.mul.eu(sigma=nysigma,X0 = 2),type="l",ylim=c(-2,2),ylab="X(t)")
plot(time,X.fun.stoch.mul.eu(sigma=nysigma,X0 = 1/2),type="l",ylim=c(-2,2),ylab="X(t)")
plot(time,X.fun.stoch.mul.eu(sigma=nysigma,X0 = -1/2),type="l",ylim=c(-2,2),ylab="X(t)")
plot(time,X.fun.stoch.mul.eu(sigma=nysigma,X0 = -2),type="l",ylim=c(-2,2),ylab="X(t)")

## Plot 4 trajectories against time for 4 different initial conditions, each
## in a different panel, simulated with the Milstein scheme:
par(mfrow=c(2,2))     ## Allows multiple panels in the same plot
plot(time,X.fun.stoch.mul.mil(sigma=nysigma,X0 = 2),type="l",ylim=c(-2,2),ylab="X(t)")
plot(time,X.fun.stoch.mul.mil(sigma=nysigma,X0 = 1/2),type="l",ylim=c(-2,2),ylab="X(t)")
plot(time,X.fun.stoch.mul.mil(sigma=nysigma,X0 = -1/2),type="l",ylim=c(-2,2),ylab="X(t)")
plot(time,X.fun.stoch.mul.mil(sigma=nysigma,X0 = -2),type="l",ylim=c(-2,2),ylab="X(t)")

some.seed <- 555

## Large step size:
dt <- 0.1

par(mfrow=c(1,1))     ## One panel in the plot
set.seed(some.seed)   ## Sets the seed of the random number generator
## Plot the trajectory, simulated with the Euler scheme:
plot(time,X.fun.stoch.mul.eu(sigma=nysigma,X0 = 2),type="l",ylim=c(-2,2),ylab="X(t)",col=2)
set.seed(some.seed)   ## Sets the seed of the random number generator
## Plot the trajectory on top, simulated with the Milstein scheme:
lines(time,X.fun.stoch.mul.mil(sigma=nysigma,X0 = 2),type="l",col=4)

## Small step size:
dt <- 0.01

par(mfrow=c(1,1))     ## One panel in the plot
set.seed(some.seed)   ## Sets the seed of the random number generator
## Plot the trajectory, simulated with the Euler scheme:
plot(time,X.fun.stoch.mul.eu(sigma=nysigma,X0 = 2),type="l",ylim=c(-2,2),ylab="X(t)",col=2)
set.seed(some.seed)   ## Sets the seed of the random number generator
## Plot the trajectory on top, simulated with the Milstein scheme:
lines(time,X.fun.stoch.mul.mil(sigma=nysigma,X0 = 2),type="l",col=4)


########################
###### EXERCISE 5 ######
########################
## NUMERICAL PROBLEMS ##
########################

## Define long time sequence:
time <- seq(0,500,dt)
N    <- length(time)

## Plot 4 trajectories against time for 4 different initial conditions:
plot(time,X.fun.det(time,X0 = 2),type="l",ylim=c(-2,2),ylab="X(t)")
lines(time,X.fun.det(time,X0 = 0.2),type="l")
lines(time,X.fun.det(time,X0 = -2),type="l")
lines(time,X.fun.det(time,X0 = -0.2),type="l")

## Problems for large t occur because the computer does not allow
## division between very large numbers. Could be solved by calculating
## the exact increments from time point to time point, instead of from
## the initial condition, such that t never becomes large in the actual
## computations (since we evaluate e^t and e^(2t))