Verifying Itô’s formula through simulation

Itô’s formula, in its simplest version, states that dYt=h(Xt)dXt+12h(Xt)(dXt)2 dY_t = h'(X_t) ~dX_t + \frac 12 h''(X_t) (dX_t)^2 when Yt=h(Xt)Y_t=h(X_t) and {Xt}\{X_t\} is an Itô process dXt=Ftdt+GtdBt. dX_t = F_t ~dt + G_t ~dB_t . Moreover, the quadratic variation [X]t[X]_t is given by d[X]t=(dXt)2=Gt2dtd[X]_t = (dX_t)^2 = G_t^2 ~dt, so dYt=h(Xt)Ftdt+h(Xt)GtdBt+12h(Xt)Gt2dt. dY_t = h'(X_t) F_t ~dt + h'(X_t) G_t ~dB_t + \frac 12 h''(X_t) G_t^2 ~dt . See (Thygesen 2023) for a description of this theoretical result.The aim of this vignette is to verify this numerically for a specific example, and to discuss discretization errors when we verify the formula in discrete time.

In our example, {Xt}\{X_t\} is stochastic logistic growth given by the Itô stochastic differential equation dXt=Xt(1Xt)dt+σXtdBt dX_t = X_t(1-X_t) ~dt + \sigma X_t ~dB_t with σ=0.25\sigma = 0.25 and the initial condition X0=x0=0.01X_0=x_0=0.01. This is known as stochastic logistic growth (Thygesen 2023). We simulate a sample path of XtX_t using the Euler-Maruyama algorithm:

## Loading required package: SDEtools
  set.seed(123456)

  f <- function(x) x*(1-x)
  g <- function(x) sigma*x
  
  t <- seq(0,10,0.01)
  B <- rBM(t)
  x0 <- 0.01
  sigma <- 0.25
  X <- euler(f,g,t,x0,B)$X
  
  plot(t,X,type="l")

We now map this process through the function h(x)=logxh(x)=\log x to generate Yt=h(Xt)Y_t=h(X_t):

  h <- function(x) log(x)
  hp <- function(x) 1/x
  hpp <- function(x) -1/x^2

  Y = h(X)
  plot(t,Y,type="l")

We now verify the first variant of Itô’s formula, viz. Yt=Y0+0th(Xs)dXs+0t12h(Xs)d[X]s. Y_t = Y_0 + \int_0^t h'(X_s) ~dX_s + \int_0^t \frac 12 h''(X_s) ~d[X]_s. We first generate the quadratic variation [X]t[X]_t, and then perform numerical Itô integration:

qvX = QuadraticVariation(X)
Y1 <- h(x0) + itointegral(hp(X),X) + 0.5*itointegral(hpp(X),qvX)

plot(t,Y,type="l")
lines(t,Y1,col="red")

The two lines lie nicely on top of each other, but there is some discretization error. We will look at this error later.

We now do the same for the second formulation, where we use the result d[X]t=g2(Xt)dtd[X]_t = g^2(X_t)~dt.

Y2 <- h(x0) + itointegral(hp(X),X) + 0.5*itointegral(hpp(X)*g(X)^2,t)

plot(t,Y,type="l")
lines(t,Y2,col="green")

Again we see that the two lines lie on top of each other, although the discretization error in this case appears to be a little more visible.

Next, we could manipulate the integral equation for {Yt}\{Y_t\} further. Inserting the forms for hh', hh'', dXtdX_t and (dXt)2(dX_t)^2, we get

dYt=(1Xt)dt+σdBt12σ2dt dY_t = (1-X_t) ~dt + \sigma ~dB_t - \frac 12 \sigma^2 ~dt

which leads to the following representation of {Yt}\{Y_t\}:

Y3 <- h(x0) + itointegral(1-X-0.5*sigma^2,t) + sigma*B

plot(t,Y,type="l")
lines(t,Y3,col="blue")

Finally, we can substitute Xt=exp(Yt)X_t = \exp(Y_t) in this expression, thus obtaining an Itô stochastic differential equation for {Yt}\{Y_t\}, viz.:

dYt=(1eYt12σ2)dt+σdBt. dY_t = (1-e^{Y_t} - \frac 12 \sigma^2) ~dt + \sigma ~dB_t.

The following code solves this equation and compares with the previous solution:

  fY <- function(y) 1-exp(y)-0.5*sigma^2
  gY <- function(y) sigma
  Y4 <- euler(fY,gY,t,log(x0),B)$X
  
  plot(t,Y,type="l")
  lines(t,Y4,col="purple")

The discretization error

We now examine the discretization error associated with the numerical Itô integration. The first formulation,

dYt=h(Xt)dXt+12h(Xt)(dXt)2 dY_t = h'(X_t) ~dX_t + \frac 12 h''(X_t) (dX_t)^2

leads to a discretization error, because we are effectively truncating at second order in Taylor’s formula when computing increments over each time step:

Yt+k=h(Xt+k)=h(Xt)+h(Xt)(Xt+kXt)+12h(Xt)(Xt+kXt)2+R Y_{t+k} = h(X_{t+k}) = h(X_t) + h'(X_t) (X_{t+k}-X_t) + \frac 12 h''(X_t) (X_{t+k}-X_t)^2 + R

so that this version of Itô’s formula ignores the higher order terms RR. We can compute the error:

  plot(t,Y-Y1,type="l")

To estimate the error, note from Taylor’s formula that the leading error term is 16h(Xt)(Xt+kXt)3\frac 16 h'''(X_t) ~(X_{t+k}-X_t)^3. Here, h(x)=2/x3h'''(x)=2/x^3. During the growth phase, the increment Xt+kXtX_{t+k}-X_t has conditional expectation kXtkX_t and variance kσ2Xt2k\sigma^2 X_t^2, so the mean cubed is Xt3(k3+3k2σ2)X_t^3(k^3+ 3 k^2\sigma^2). As a result, the error should grow with rate (k3+3k2σ2)/36.6106(k^3+3k^2\sigma^2)/3 \approx 6.6\cdot 10^{-6} per time step, i.e. with a rate 6.61046.6\cdot 10^{-4}. So after 5 time units, the error should be 0.00330.0033. We see that there is good agreement with the simulation result.

We now consider the discretization error in the second formulation, dYt=h(Xt)dXt+12h(Xt)g2(Xt)dt. dY_t = h'(X_t) ~dX_t + \frac 12 h''(X_t) g^2(X_t) ~dt . Here, there is also the effect that the increment in discretized quadratic variation [X]t[X]_t, i.e. (Xt+kXt)2(X_{t+k}-X_t)^2, does not exactly equal g2(Xt)kg^2(X_t) ~k. In fact, the conditional variance of the increment equals g2(Xt)kg^2(X_t) ~k, but in discrete time, there is both a random component and a mean. This error can be traced back the underlying Brownian motion: In discrete time, its quadratic variation does not exactly equal time. In this case, this error turns out to be larger than the one we examined previously:

  plot(t,Y-Y2,type="l")

To assess this error coarsely, note that it comes from the term 12h(Xt)[(Xt+kXt)2g2(Xt)k] \frac 12 h''(X_t) [ (X_{t+k}-X_t)^2 - g^2(X_t) ~k] and the conditional expectation of this term is k2/2-k^2/2 per time step during the growth phase, i.e. it grows with rate k/2-k/2. So after 5 time units, we expect that this error has grown to 0.025. We see a good agreement with the simulation, even if the variation between realizations are greater.

This vignette has hopefully illuminated both where Itô’s formula comes from, and the errors that may arise when we investigate stochastic differential equation in discrete time.

References

Thygesen, Uffe Høgsbro. 2023. Stochastic Differential Equations for Science and Engineering. Taylor & Francis. http://uffe-h-thygesen.github.io.