--- title: "Parametric Estimation of 1-D Stochastic Differential Equation" author: - A.C. Guidoum^[Department of Mathematics and Computer Science, Faculty of Sciences and Technology, University of Tamanghasset, Algeria, E-mail (acguidoum@univ-tam.dz)] and K. Boukhetala^[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)] date: "`r Sys.Date()`" output: knitr:::html_vignette: toc: yes vignette: > %\VignetteIndexEntry{Parametric Estimation of 1-D Stochastic Differential Equation} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, echo = F, message = F, results = 'hide',screenshot.force=FALSE} library(Sim.DiffProc) library(knitr) knitr::opts_chunk$set(comment="",prompt=TRUE, fig.show='hold', warning=FALSE, message=FALSE) options(prompt="R> ",scipen=16,digits=5,warning=FALSE, message=FALSE, width = 70) ``` # The `fitsde()` function The [Sim.DiffProc](https://CRAN.R-project.org/package=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: - `data` a univariate time series (`ts` object). - `drift` and `diffusion` indicate drift and diffusion coefficient of the SDE, is an `expression` of two variables `t`, `x` and `theta` names of the parameters, and must be nominated by a vector of `theta = (theta[1], theta[2],..., theta[p])` for reasons of symbolic derived in approximation methods. - `start` must be specified as a named list, where the names of the elements of the list correspond to the names of the parameters as they appear in the `drift` and `diffusion` coefficient. - The `pmle` argument must be a `character` string specifying the method to use, can be either: `"euler"` Euler pseudo-likelihood, `"ozaki"` Ozaki pseudo-likelihood, `"shoji"` Shoji pseudo-likelihood and `"kessler"` Kessler pseudo-likelihood. - `optim.method` select the optimization method (`"L-BFGS-B"` is used by default), and further arguments to pass to `optim` function. - `lower` and `upper` bounds on the variables for the `Brent` or `L-BFGS-B` method. 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: - `coef`: which extracts model coefficients from objects returned by `fitsde`. - `vcov`: returns the variance-covariance matrix of the parameters of a fitted model objects. - `logLik`: extract log-likelihood. - `AIC`: calculating Akaike's Information Criterion for fitted model objects. - `BIC`: calculating Schwarz's Bayesian Criterion for fitted model objects. - `confint`: computes confidence intervals for one or more parameters in a fitted model objects. 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} set.seed(12345, kind = "L'Ecuyer-CMRG") f <- expression( (1+2*x) ) ; g <- expression( 0.5*x^0.3 ) sim <- snssde1d(drift=f,diffusion=g,x0=2,N=10^4,Dt=10^-4) 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, message=FALSE, warning=FALSE} fx <- expression( theta[1]+theta[2]*x ) ## drift coefficient of model gx <- expression( theta[3]*x^theta[4] ) ## diffusion coefficient of model fitmod <- fitsde(data = mydata, drift = fx, diffusion = gx, start = list(theta1=1, theta2=1, theta3=1,theta4=1),pmle="euler") fitmod ``` The estimated coefficients are extracted from the output object `fitmod` as follows: ```{r} coef(fitmod) ``` We can use the `summary` function to produce result summaries of output object: ```{r} summary(fitmod) ``` `vcov` for variance-covariance matrice, and extract log-likelihood by `logLik`: ```{r} vcov(fitmod) logLik(fitmod) AIC(fitmod) BIC(fitmod) ``` Computes confidence intervals for one or more parameters in a fitted SDE: ```{r} confint(fitmod, level=0.95) ``` ## 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} set.seed(1234, kind = "L'Ecuyer-CMRG") f <- expression( 3*(2-x) ) ; g <- expression( 0.5 ) sim <- snssde1d(drift=f,diffusion=g,x0=5,Dt=0.01) 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, message=FALSE, warning=FALSE} fx <- expression( theta[1]*(theta[2]- x) ) ## drift coefficient of model gx <- expression( theta[3] ) ## diffusion coefficient of model fitmod <- fitsde(data=HWV,drift=fx,diffusion=gx,start = list(theta1=1,theta2=1, theta3=1),pmle="ozaki") summary(fitmod) ``` 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) ``` ## 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} set.seed(1234, kind = "L'Ecuyer-CMRG") f <- expression(-2*x*t) ; g <- expression(0.2*x) sim <- snssde1d(drift=f,diffusion=g,N=1000,Dt=0.001,x0=10) 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, message=FALSE, warning=FALSE} fx <- expression( theta[1]*x*t ) ## drift coefficient of model gx <- expression( theta[2]*x ) ## diffusion coefficient of model fitmod <- fitsde(data=mydata,drift=fx,diffusion=gx,start = list(theta1=1, theta2=1),pmle="shoji",lower=c(-3,0),upper=c(-1,1)) summary(fitmod) ``` ## 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} set.seed(1234, kind = "L'Ecuyer-CMRG") f <- expression(3*t*(sqrt(t)-x)) ; g <- expression(0.3*t) sim <- snssde1d(drift=f,diffusion=g,M=1,N=1000,x0=2,Dt=0.001) 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, message=FALSE, warning=FALSE} fx <- expression( theta[1]*t* ( theta[2]*sqrt(t) - x ) ) ## drift coefficient of model gx <- expression( theta[3]*t ) ## diffusion coefficient of model fitmod <- fitsde(data=mydata,drift=fx,diffusion=gx,start = list(theta1=1, theta2=1,theta3=1),pmle="kessler") summary(fitmod) ``` # 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} set.seed(1234, kind = "L'Ecuyer-CMRG") f <- expression( 2*x ) g <- expression( 0.3*x^0.5 ) sim <- snssde1d(drift=f,diffusion=g,M=1,N=10000,x0=2,Dt=0.0001) mydata <- sim$X ``` We test the performance of the AIC statistics for the four competing models, we can proceed as follows: ```{r, message=FALSE, warning=FALSE} ## True model fx <- expression( theta[1]*x ) gx <- expression( theta[2]*x^theta[3] ) truemod <- fitsde(data=mydata,drift=fx,diffusion=gx,start = list(theta1=1, theta2=1,theta3=1),pmle="euler") ## competing model 1 fx1 <- expression( theta[1]+theta[2]*x ) gx1 <- expression( theta[3]*x^theta[4] ) mod1 <- fitsde(data=mydata,drift=fx1,diffusion=gx1,start = list(theta1=1, theta2=1,theta3=1,theta4=1),pmle="euler") ## competing model 2 fx2 <- expression( theta[1]+theta[2]*x ) gx2 <- expression( theta[3]*sqrt(x) ) mod2 <- fitsde(data=mydata,drift=fx2,diffusion=gx2,start = list(theta1=1, theta2=1,theta3=1),pmle="euler") ## competing model 3 fx3 <- expression( theta[1] ) gx3 <- expression( theta[2]*x^theta[3] ) mod3 <- fitsde(data=mydata,drift=fx3,diffusion=gx3,start = list(theta1=1, theta2=1,theta3=1),pmle="euler") ## Computes AIC AIC <- c(AIC(truemod),AIC(mod1),AIC(mod2),AIC(mod3)) Test <- data.frame(AIC,row.names = c("True mod","Comp mod1","Comp mod2","Comp mod3")) Bestmod <- rownames(Test)[which.min(Test[,1])] Bestmod ``` the estimates under the different models: ```{r} Theta1 <- c(coef(truemod)[[1]],coef(mod1)[[1]],coef(mod2)[[1]],coef(mod3)[[1]]) Theta2 <- c(coef(truemod)[[2]],coef(mod1)[[2]],coef(mod2)[[2]],coef(mod3)[[2]]) Theta3 <- c(coef(truemod)[[3]],coef(mod1)[[3]],coef(mod2)[[3]],coef(mod3)[[3]]) Theta4 <- c("",round(coef(mod1)[[4]],5),"","") Parms <- data.frame(Theta1,Theta2,Theta3,Theta4,row.names = c("True mod", "Comp mod1","Comp mod2","Comp mod3")) Parms ``` ## 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](https://cran.r-project.org/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](https://cran.r-project.org/package=yuima) package, here we confirm the results of the estimates by several approximation methods. ```{r 01,fig.env='figure*', fig.cap=' The U.S. Interest Rates monthly form $06/1964$ to $12/1989$ ',fig.width=6,fig.height=4} data(Irates) rates <- Irates[, "r1"] X <- window(rates, start = 1964.471, end = 1989.333) plot(X) ``` we can now use all previous methods by `fitsde` function to estimate the parameters of CKLS model as follows: ```{r, message=FALSE, warning=FALSE} fx <- expression( theta[1]+theta[2]*x ) ## drift coefficient of CKLS model gx <- expression( theta[3]*x^theta[4] ) ## diffusion coefficient of CKLS model pmle <- eval(formals(fitsde.default)$pmle) 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))) Coef <- data.frame(do.call("cbind",lapply(1:4,function(i) coef(fitres[[i]])))) 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) names(Coef) <- c(pmle) names(Info) <- c("logLik","AIC","BIC") Coef Info ``` 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 02,fig.env='figure*', fig.cap='The path mean of the solution of the CKLS model with the estimated parameters and real data ',fig.width=6,fig.height=4} set.seed(1234, kind = "L'Ecuyer-CMRG") f <- expression( (2.076-0.263*x) ) g <- expression( 0.130*x^1.451 ) mod <- snssde1d(drift=f,diffusion=g,x0=X[1],M=500, N=length(X),t0=1964.471, T=1989.333) mod plot(mod,type="n",ylim=c(0,30)) lines(X,col=4,lwd=2) lines(time(mod),apply(mod$X,1,mean),col=2,lwd=2) lines(time(mod),apply(mod$X,1,bconfint,level=0.95)[1,],col=5,lwd=2) lines(time(mod),apply(mod$X,1,bconfint,level=0.95)[2,],col=5,lwd=2) 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](snssde.html). 2. [`bridgesdekd()` & `dsdekd()` & `rsdekd()` - Constructs and Analysis of Bridges Stochastic Differential Equations](bridgesde.html). 3. [`fptsdekd()` & `dfptsdekd()` - Monte-Carlo Simulation and Kernel Density Estimation of First passage time](fptsde.html). 4. [`MCM.sde()` & `MEM.sde()` - Parallel Monte-Carlo and Moment Equations for SDEs](mcmsde.html). 5. [`TEX.sde()` - Converting Sim.DiffProc Objects to LaTeX](sdetotex.html). 6. [`fitsde()` - Parametric Estimation of 1-D Stochastic Differential Equation](fitsde.html). # 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 https://www.jstatsoft.org/v57/i04. 2. Guidoum AC, Boukhetala K (2020). "Performing Parallel Monte Carlo and Moment Equations Methods for Itô and Stratonovich Stochastic Differential Systems: R Package Sim.DiffProc". Journal of Statistical Software, 96(2), 1--82. https://doi.org/10.18637/jss.v096.i02 3. Iacus SM (2008). Simulation and Inference for Stochastic Differential Equations: With R Examples. Springer-Verlag, New York.