Parametric Estimation of 1-D Stochastic Differential Equation

A.C. Guidoum1 and K. Boukhetala2

2018-10-17

The fitsde() function

The Sim.DiffProc package implements pseudo-maximum likelihood via the fitsde() function. The interface and the output of the fitsde() function are made as similar as possible to those of the standard mle function in the stats4 package of the basic R system. The main arguments to fitsde consist:

The functions of type S3 method (as similar of the standard mle function in the stats4 package of the basic R system for the class fitsde are the following:

The following we explain how to use this function to estimate a SDE with different approximation methods.

Euler method

Consider a process solution of the general stochastic differential equation: \[\begin{equation}\label{eq02} dX_{t}= f(X_{t},\underline{\theta}) dt + g(X_{t},\underline{\theta}) dW_{t}, \quad \quad t \geq 0 \, , X_{0} = x_{0}, \end{equation}\] The Euler scheme produces the discretization (\(\Delta t \rightarrow 0\)): \[\begin{equation*} X_{t+\Delta t} - X_{t} = f(X_{t},\theta) \Delta t+ g(X_{t},\theta) (W_{t+\Delta t} -W_{t}), \end{equation*}\] The increments \(X_{t+\Delta t} - X_{t}\) are then independent Gaussian random variables with mean: \(\text{E}_{x} = f(X_{t},\theta)\Delta t\), and variance: \(\text{V}_{x} = g^{2}(X_{t},\theta) \Delta t\). Therefore the transition density of the process can be written as: \[\begin{equation*} p_{\theta}(t,y|x)=\frac{1}{\sqrt{2\pi t g^{2}(x,\theta)}} \exp\left(-\frac{\left(y-x-f(x,\theta)t\right)^2}{2tg^{2}(x,\theta)}\right) \end{equation*}\] Then the log-likelihood is: \[\begin{equation}\label{eq08} h_{n}(\theta|X_{1},X_{2},\dots,X_{n})=-\frac{1}{2}\left(\sum_{i=1}^{n} \frac{(X_{i}-X_{i-1}-f(X_{i-1},\theta)\Delta)^2}{\sigma^2 \Delta t} + n \log(2\pi \sigma^2 \Delta t)\right) \end{equation}\] As an example, we consider the Chan-Karolyi-Longstaff-Sanders (CKLS) model: \[\begin{equation}\label{eq09} dX_{t} = (\theta_{1}+\theta_{2} X_{t}) dt + \theta_{3} X_{t}^{\theta_{4}} dW_{t},\qquad X_{0}=2 \end{equation}\]

with \(\theta_{1}=1\), \(\theta_{2}=2\), \(\theta_{3}=0.5\) and \(\theta_{4}=0.3\). Before calling fitsde, we generate sampled data \(X_{t_{i}}\), with \(\Delta t =10^{-4}\), as following:

R> f <- expression( (1+2*x) ) ; g <- expression( 0.5*x^0.3 )
R> sim    <- snssde1d(drift=f,diffusion=g,x0=2,N=10^4,Dt=10^-4)
R> mydata <- sim$X

we set the initial values for the optimizer as \(\theta_{1}=\theta_{2}=\theta_{3}=\theta_{4}=1\), and we specify the coefficients drift and diffusion as expressions. you can use the upper and lower limits of the search region used by the optimizer (here using the default method "L-BFGS-B"), and specifying the method to use with pmle="euler".

R> fx <- expression( theta[1]+theta[2]*x ) ## drift coefficient of model
R> gx <- expression( theta[3]*x^theta[4] ) ## diffusion coefficient of model 
R> fitmod <- fitsde(data = mydata, drift = fx, diffusion = gx, start = list(theta1=1, theta2=1,
+                 theta3=1,theta4=1),pmle="euler")
R> fitmod

Call:
fitsde(data = mydata, drift = fx, diffusion = gx, start = list(theta1 = 1, 
    theta2 = 1, theta3 = 1, theta4 = 1), pmle = "euler")

Coefficients:
 theta1  theta2  theta3  theta4 
1.11883 2.08943 0.50831 0.29491 

The estimated coefficients are extracted from the output object fitmod as follows:

R> coef(fitmod)
 theta1  theta2  theta3  theta4 
1.11883 2.08943 0.50831 0.29491 

We can use the summary function to produce result summaries of output object:

R> summary(fitmod)
Pseudo maximum likelihood estimation

Method:  Euler
Call:
fitsde(data = mydata, drift = fx, diffusion = gx, start = list(theta1 = 1, 
    theta2 = 1, theta3 = 1, theta4 = 1), pmle = "euler")

Coefficients:
       Estimate Std. Error
theta1  1.11883   1.493710
theta2  2.08943   0.198762
theta3  0.50831   0.010823
theta4  0.29491   0.010796

-2 log L: -66286 

vcov for variance-covariance matrice, and extract log-likelihood by logLik:

R> vcov(fitmod)
             theta1        theta2        theta3        theta4
theta1  2.231169308 -0.2443954942 -0.0000339028  0.0000462854
theta2 -0.244395494  0.0395063313  0.0000067968 -0.0000083309
theta3 -0.000033903  0.0000067968  0.0001171369 -0.0001102133
theta4  0.000046285 -0.0000083309 -0.0001102133  0.0001165533
R> logLik(fitmod)
[1] 33143
R> AIC(fitmod)
[1] -66278
R> BIC(fitmod)
[1] -66267

Computes confidence intervals for one or more parameters in a fitted SDE:

R> confint(fitmod, level=0.95)
          2.5 %  97.5 %
theta1 -1.80879 4.04645
theta2  1.69986 2.47900
theta3  0.48709 0.52952
theta4  0.27375 0.31607

Ozaki method

The second approach we present is the Ozaki method, and it works for homogeneous stochastic differential equations. Consider the stochastic differential equation: \[\begin{equation}\label{eq10} dX_{t}= f(X_{t},\underline{\theta}) dt + \sigma dW_{t}, \quad \quad t \geq 0 \, , X_{0} = x_{0}, \end{equation}\] where \(\sigma\) is supposed to be constant. We just recall that the transition density for the Ozaki method is Gaussian, we have that: \(X_{t+\Delta t}|X_{t} = x \sim \mathcal{N}(\text{E}_{x},\text{V}_{x})\), where: \[\begin{align}\label{eq11} \text{E}_{x} &= x + \frac{f(x)}{\partial_{x} f(x)} \left( e^{\partial_{x} f(x) \Delta t} - 1 \right), \quad\text{and}\quad \text{V}_{x} &= \sigma^{2} \frac{e^{2K_{x} \Delta t} -1}{2K_{x}}, \end{align}\] with: \[\begin{equation*} K_{x} = \frac{1}{\Delta t} \log \left(1+\frac{f(x)}{x\partial_{x}f(x)}\left(e^{\partial_{x}f(x) \Delta t}-1\right) \right) \end{equation*}\]

It is always possible to transform process \(X_t\) with a constant diffusion coefficient using the Lamperti transform.

Now we consider the Vasicek model, with \(f(x,\theta) = \theta_{1} (\theta_{2}- x)\) and constant volatility \(g(x,\theta) = \theta_{3}\), \[\begin{equation}\label{eq12} dX_{t} = \theta_{1} (\theta_{2}- X_{t}) dt + \theta_{3} dW_{t},\qquad X_{0}=5 \end{equation}\]

with \(\theta_{1}=3\), \(\theta_{2}=2\) and \(\theta_{3}=0.5\), we generate sampled data \(X_{t_{i}}\), with \(\Delta t =10^{-2}\), as following:

R> f <- expression( 3*(2-x) ) ; g <- expression( 0.5 )
R> sim <- snssde1d(drift=f,diffusion=g,x0=5,Dt=0.01)
R> HWV <- sim$X

we set the initial values for the optimizer as \(\theta_{1}=\theta_{2}=\theta_{3}=1\), and we specify the coefficients drift and diffusion as expressions. Specifying the method to use with pmle="ozaki", which can easily be implemented in R as follows:

R> fx <- expression( theta[1]*(theta[2]- x) ) ## drift coefficient of model 
R> gx <- expression( theta[3] )           ## diffusion coefficient of model 
R> fitmod <- fitsde(data=HWV,drift=fx,diffusion=gx,start = list(theta1=1,theta2=1,
+                   theta3=1),pmle="ozaki")
R> summary(fitmod)
Pseudo maximum likelihood estimation

Method:  Ozaki
Call:
fitsde(data = HWV, drift = fx, diffusion = gx, start = list(theta1 = 1, 
    theta2 = 1, theta3 = 1), pmle = "ozaki")

Coefficients:
       Estimate Std. Error
theta1  3.47753   0.405022
theta2  1.97249   0.048249
theta3  0.50968   0.011398

-2 log L: -3115.9 

If you want to have confidence intervals \(\theta_{1}\) and \(\theta_{2}\) parameters only, using the command confint as following:

R> confint(fitmod,parm=c("theta1","theta2"),level=0.95)
        2.5 % 97.5 %
theta1 2.6837 4.2714
theta2 1.8779 2.0671

Shoji-Ozaki method

An extension of the method to Ozaki the more general case where the drift is allowed to depend on the time variable \(t\), and also the diffusion coefficient can be varied is the method Shoji and Ozaki. Consider the stochastic differential equation: \[\begin{equation}\label{eq13} dX_{t}= f(t,X_{t},\underline{\theta}) dt + g(X_{t},\underline{\theta}) dW_{t}, \quad \quad t \geq 0 \, , X_{0} = x_{0}, \end{equation}\] the transition density for the Shoji-Ozaki method is Gaussian, we have that: \(X_{t+\Delta t}|X_{t} = x \sim \mathcal{N}\left(\mathrm{A}_{(t,x)}x,\mathrm{B}^{2}_{(t,x)}\right)\), where: \[\begin{align}\label{eq14} \mathrm{A}_{(t,x)} &= 1+ \frac{f(t,x)}{x\mathrm{L}_{t}} \left(e^{\mathrm{L}_{t}\Delta t }-1\right)+\frac{\mathrm{M}_{t}}{x\mathrm{L}^{2}_{t}} \left(e^{\mathrm{L}_{t} \Delta t}-1-\mathrm{L}_{t}\Delta t\right), \\ \mathrm{B}_{(t,x)} &= g(x) \sqrt{\frac{e^{2\mathrm{L}_{t} \Delta t}-1}{2\mathrm{L}_{t}}}, \end{align}\] with: \[\begin{equation*} \mathrm{L}_{t} = \partial_{x} f(t,x) \quad\text{and}\quad \mathrm{M}_{t} = \frac{g^{2}(x)}{2} \partial_{xx} f(t,x)+ \partial_{t} f(t,x). \end{equation*}\] As an example, we consider the following model: \[\begin{equation}\label{eq15} dX_{t} = a(t)X_{t} dt + \theta_{2}X_{t} dW_{t},\qquad X_{0}=10 \end{equation}\]

with: \(a(t) = \theta_{1}t\), and we generate sampled data \(X_{t_{i}}\), with \(\theta_{1}=-2\), \(\theta_{2}=0.2\) and time step \(\Delta t =10^{-3}\), as following:

R> f <- expression(-2*x*t) ; g <- expression(0.2*x)
R> sim <- snssde1d(drift=f,diffusion=g,N=1000,Dt=0.001,x0=10)
R> mydata <- sim$X

we set the initial values for the optimizer as \(\theta_{1}=\theta_{2}=1\), and we specify the method to use with pmle="shoji":

R> fx <- expression( theta[1]*x*t ) ## drift coefficient of model 
R> gx <- expression( theta[2]*x )   ## diffusion coefficient of model 
R> fitmod <- fitsde(data=mydata,drift=fx,diffusion=gx,start = list(theta1=1,
+                  theta2=1),pmle="shoji",lower=c(-3,0),upper=c(-1,1))
R> summary(fitmod)
Pseudo maximum likelihood estimation

Method:  Shoji
Call:
fitsde(data = mydata, drift = fx, diffusion = gx, start = list(theta1 = 1, 
    theta2 = 1), pmle = "shoji", lower = c(-3, 0), upper = c(-1, 
    1))

Coefficients:
       Estimate Std. Error
theta1 -2.66079  0.3483607
theta2  0.20094  0.0044949

-2 log L: -3500.5 

Kessler method

Kessler (1997) proposed to use a higher-order Ito-Taylor expansion to approximate the mean and variance in a conditional Gaussian density. Consider the stochastic differential equation \(dX_{t}= f(X_{t},\underline{\theta}) dt + g(X_{t},\underline{\theta}) dW_{t}\) , the transition density by Kessler method is: \(X_{t+\Delta t}|X_{t} = x \sim \mathcal{N}\left(\text{E}_{x},\text{V}_{x}\right)\), where: \[\begin{align}\label{eq16} \text{E}_{x} &= x + f(t,x) \Delta t+\left(f(t,x)\partial_{x}f(t,x) + \frac{1}{2} g^{2}(t,x) \partial_{xx}g(t,x)\right)\frac{(\Delta t)^2}{2}, \\ \text{V}_{x} &= x^2 +(2f(t,x)x+g^{2}(t,x)) \Delta t +\bigg(2f(t,x)\left(\partial_{x}f(t,x)x+f(t,x)+g(t,x)\partial_{x}g(t,x)\right) \nonumber\\ &\quad+g^{2}(t,x)\left(\partial_{xx}f(t,x)x+2\partial_{x}f(t,x)+\partial_{x}g^{2}(t,x)+g(t,x)\partial_{xx}g(t,x)\right)\bigg)\frac{(\Delta t)^2}{2}-\text{E}^{2}_{x}. \end{align}\] We consider the following Hull-White (extended Vasicek) model: \[\begin{equation}\label{eq17} dX_{t} = a(t)(b(t)-X_{t}) dt + \sigma(t) dW_{t},\qquad X_{0}=2 \end{equation}\]

with: \(a(t) = \theta_{1}t\) and \(b(t)=\theta_{2}\sqrt{t}\), the volatility depends on time: \(\sigma(t)=\theta_{3}t\). We generate sampled data of \(X_t\), with \(\theta_{1}=3\), \(\theta_{2}=1\) and \(\theta_{3}=0.3\), time step \(\Delta t =10^{-3}\), as following:

R> f <- expression(3*t*(sqrt(t)-x)) ; g <- expression(0.3*t)
R> sim <- snssde1d(drift=f,diffusion=g,M=1,N=1000,x0=2,Dt=0.001)
R> mydata <- sim$X

we set the initial values for the optimizer as \(\theta_{1}=\theta_{2}=\theta_{3}=1\), and we specify the method to use with pmle="kessler":

R> fx <- expression( theta[1]*t* ( theta[2]*sqrt(t) - x ) ) ## drift coefficient of model 
R> gx <- expression( theta[3]*t ) ## diffusion coefficient of model 
R> fitmod <- fitsde(data=mydata,drift=fx,diffusion=gx,start = list(theta1=1,
+                   theta2=1,theta3=1),pmle="kessler")
R> summary(fitmod)
Pseudo maximum likelihood estimation

Method:  Kessler
Call:
fitsde(data = mydata, drift = fx, diffusion = gx, start = list(theta1 = 1, 
    theta2 = 1, theta3 = 1), pmle = "kessler")

Coefficients:
       Estimate Std. Error
theta1  3.06677  0.3854323
theta2  1.18670  0.1686620
theta3  0.29689  0.0066424

-2 log L: -8484.9 

The fitsde() in practice

Model selection via AIC

Let the following models: \[\begin{align*} % \nonumber to remove numbering (before each equation) dX_{t} &= \theta_{1} X_{t} dt + \theta_{2} X_{t}^{\theta_{3}} dW_{t}, &\text{(true model)}\\ dX_{t} &= (\theta_{1}+\theta_{2} X_{t}) dt + \theta_{3} X_{t}^{\theta_{4}} dW_{t},&\text{(competing model 1)}\\ dX_{t} &= (\theta_{1}+\theta_{2} X_{t}) dt + \theta_{3} \sqrt{X_{t}} dW_{t}, &\text{(competing model 2)}\\ dX_{t} &= \theta_{1} dt + \theta_{2} X_{t}^{\theta_{3}} dW_{t}, &\text{(competing model 3)} \end{align*}\]

We generate data from true model with parameters \(\underline{\theta}=(2,0.3,0.5)\), initial value \(X_{0}=2\) and \(\Delta t =10^{-4}\), as following:

R> f <- expression( 2*x )
R> g <- expression( 0.3*x^0.5 )
R> sim <- snssde1d(drift=f,diffusion=g,M=1,N=10000,x0=2,Dt=0.0001)
R> mydata <- sim$X

We test the performance of the AIC statistics for the four competing models, we can proceed as follows:

R> ## True model
R> fx <- expression( theta[1]*x )
R> gx <- expression( theta[2]*x^theta[3] )
R> truemod <- fitsde(data=mydata,drift=fx,diffusion=gx,start = list(theta1=1,
+                    theta2=1,theta3=1),pmle="euler")
R> ## competing model 1
R> fx1 <- expression( theta[1]+theta[2]*x )
R> gx1 <- expression( theta[3]*x^theta[4] )
R> mod1 <- fitsde(data=mydata,drift=fx1,diffusion=gx1,start = list(theta1=1,
+           theta2=1,theta3=1,theta4=1),pmle="euler")
R> ## competing model 2
R> fx2 <- expression( theta[1]+theta[2]*x )
R> gx2 <- expression( theta[3]*sqrt(x) )
R> mod2 <- fitsde(data=mydata,drift=fx2,diffusion=gx2,start = list(theta1=1,
+            theta2=1,theta3=1),pmle="euler")
R> ## competing model 3
R> fx3 <- expression( theta[1] )
R> gx3 <- expression( theta[2]*x^theta[3] )
R> mod3 <- fitsde(data=mydata,drift=fx3,diffusion=gx3,start = list(theta1=1,
+            theta2=1,theta3=1),pmle="euler")
R> ## Computes AIC
R> AIC <- c(AIC(truemod),AIC(mod1),AIC(mod2),AIC(mod3))
R> Test <- data.frame(AIC,row.names = c("True mod","Comp mod1","Comp mod2","Comp mod3"))
R> Bestmod <- rownames(Test)[which.min(Test[,1])]
R> Bestmod
[1] "Comp mod2"

the estimates under the different models:

R> Theta1 <- c(coef(truemod)[[1]],coef(mod1)[[1]],coef(mod2)[[1]],coef(mod3)[[1]])
R> Theta2 <- c(coef(truemod)[[2]],coef(mod1)[[2]],coef(mod2)[[2]],coef(mod3)[[2]])
R> Theta3 <- c(coef(truemod)[[3]],coef(mod1)[[3]],coef(mod2)[[3]],coef(mod3)[[3]])
R> Theta4 <- c("",round(coef(mod1)[[4]],5),"","")
R> Parms  <- data.frame(Theta1,Theta2,Theta3,Theta4,row.names = c("True mod",
+                       "Comp mod1","Comp mod2","Comp mod3"))
R> Parms
          Theta1  Theta2  Theta3  Theta4
True mod  1.8893 0.30760 0.48356        
Comp mod1 1.1148 1.72367 0.30745 0.48379
Comp mod2 1.9295 1.60008 0.29889        
Comp mod3 9.8694 0.30658 0.48681        

Application to real data

We make use of real data of the U.S. Interest Rates monthly form \(06/1964\) to \(12/1989\) (see Figure 1) available in package Ecdat, and we estimate the parameters \(\underline{\theta}=(\theta_{1},\theta_{2},\theta_{3},\theta_{4})\) of CKLS model. These data have been analyzed by Brouste et all (2014) with yuima package, here we confirm the results of the estimates by several approximation methods.

R> data(Irates)
R> rates <- Irates[, "r1"]
R> X <- window(rates, start = 1964.471, end = 1989.333)
R> plot(X)

we can now use all previous methods by fitsde function to estimate the parameters of CKLS model as follows:

R> fx <- expression( theta[1]+theta[2]*x ) ## drift coefficient of CKLS model
R> gx <- expression( theta[3]*x^theta[4] ) ## diffusion coefficient of CKLS model
R> pmle <- eval(formals(fitsde.default)$pmle)
R> fitres <- lapply(1:4, function(i) fitsde(X,drift=fx,diffusion=gx,pmle=pmle[i],
+                   start = list(theta1=1,theta2=1,theta3=1,theta4=1)))
R> Coef <- data.frame(do.call("cbind",lapply(1:4,function(i) coef(fitres[[i]]))))
R> Info <- data.frame(do.call("rbind",lapply(1:4,function(i) logLik(fitres[[i]]))),
+                    do.call("rbind",lapply(1:4,function(i) AIC(fitres[[i]]))),
+                    do.call("rbind",lapply(1:4,function(i) BIC(fitres[[i]]))),
+                    row.names=pmle)
R> names(Coef) <- c(pmle)
R> names(Info) <- c("logLik","AIC","BIC")
R> Coef
          euler  kessler    ozaki    shoji
theta1  2.07695  2.14335  2.11532  2.10150
theta2 -0.26319 -0.27434 -0.26905 -0.26647
theta3  0.13022  0.12598  0.12652  0.13167
theta4  1.45132  1.46917  1.46491  1.45131
R> Info
         logLik    AIC    BIC
euler   -237.88 483.76 487.15
kessler -237.78 483.57 486.96
ozaki   -237.84 483.67 487.07
shoji   -237.88 483.76 487.15

In Figure 2 we reports the sample mean of the solution of the CKLS model with the estimated parameters and real data, their empirical \(95\%\) confidence bands (from the \(2.5th\) to the \(97.5th\) percentile), we can proceed as follows:

R> f <- expression( (2.076-0.263*x) )
R> g <- expression( 0.130*x^1.451 )
R> mod <- snssde1d(drift=f,diffusion=g,x0=X[1],M=500, N=length(X),t0=1964.471, T=1989.333)
R> mod
Itô Sde 1D:
    | dX(t) = (2.076 - 0.263 * X(t)) * dt + 0.13 * X(t)^1.451 * dW(t)
Method:
    | Euler scheme with order 0.5
Summary:
    | Size of process   | N  = 299.
    | Number of simulation  | M  = 500.
    | Initial value     | x0 = 3.317.
    | Time of process   | t in [1964.5,1989.3].
    | Discretization    | Dt = 0.08343.
R> plot(mod,type="n",ylim=c(0,30))
R> lines(X,col=4,lwd=2)
R> lines(time(mod),apply(mod$X,1,mean),col=2,lwd=2)
R> lines(time(mod),apply(mod$X,1,bconfint,level=0.95)[1,],col=5,lwd=2)
R> lines(time(mod),apply(mod$X,1,bconfint,level=0.95)[2,],col=5,lwd=2)
R> legend("topleft",c("real data","mean path",paste("bound of", 95,"% confidence")),inset = .01,col=c(4,2,5),lwd=2,cex=0.8)

Further reading

  1. snssdekd() & dsdekd() & rsdekd()- Monte-Carlo Simulation and Analysis of Stochastic Differential Equations.
  2. bridgesdekd() & dsdekd() & rsdekd() - Constructs and Analysis of Bridges Stochastic Differential Equations.
  3. fptsdekd() & dfptsdekd() - Monte-Carlo Simulation and Kernel Density Estimation of First passage time.
  4. MCM.sde() & MEM.sde() - Parallel Monte-Carlo and Moment Equations for SDEs.
  5. TEX.sde() - Converting Sim.DiffProc Objects to LaTeX.
  6. fitsde() - Parametric Estimation of 1-D Stochastic Differential Equation.

References

  1. Brouste A, Fukasawa M, Hino H, Iacus SM, Kamatani K, Koike Y, Masuda H, Nomura R,Ogihara T, Shimuzu Y, Uchida M, Yoshida N (2014). The YUIMA Project: A ComputationalFramework for Simulation and Inference of Stochastic Differential Equations." Journal of Statistical Software, 57(4), 1-51. URL http://www.jstatsoft.org/v57/i04.

  2. Guidoum AC, Boukhetala K (2018). Sim.DiffProc: Simulation of Diffusion Processes. R package version 4.2, URL https://cran.r-project.org/package=Sim.DiffProc.

  3. Iacus SM (2008). Simulation and Inference for Stochastic Differential Equations: With R Examples. Springer-Verlag, New York.


  1. Department of Probabilities & Statistics, Faculty of Mathematics, University of Science and Technology Houari Boumediene, BP 32 El-Alia, U.S.T.H.B, Algeria, E-mail (acguidoum@usthb.dz)

  2. Faculty of Mathematics, University of Science and Technology Houari Boumediene, BP 32 El-Alia, U.S.T.H.B, Algeria, E-mail (kboukhetala@usthb.dz)