Nothing
Parametric Estimation of 1-D Stochastic Differential Equation"
In Sim.DiffProc: Simulation of Diffusion Processes
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 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:
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".
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:
coef(fitmod)
We can use the summary function to produce result summaries of output object:
summary(fitmod)
vcov for variance-covariance matrice, and extract log-likelihood by logLik:
vcov(fitmod)
logLik(fitmod)
AIC(fitmod)
BIC(fitmod)
Computes confidence intervals for one or more parameters in a fitted SDE:
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:
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:
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:
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:
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":
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:
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":
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:
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:
## 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:
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, 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.
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:
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:
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
snssdekd() & dsdekd() & rsdekd()- Monte-Carlo Simulation and Analysis of Stochastic Differential Equations.bridgesdekd() & dsdekd() & rsdekd() - Constructs and Analysis of Bridges Stochastic Differential Equations.fptsdekd() & dfptsdekd() - Monte-Carlo Simulation and Kernel Density Estimation of First passage time.MCM.sde() & MEM.sde() - Parallel Monte-Carlo and Moment Equations for SDEs.TEX.sde() - Converting Sim.DiffProc Objects to LaTeX.fitsde() - Parametric Estimation of 1-D Stochastic Differential Equation.
References
-
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.
-
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
-
Iacus SM (2008). Simulation and Inference for Stochastic Differential Equations: With R Examples. Springer-Verlag, New York.
Try the Sim.DiffProc package in your browser
Any scripts or data that you put into this service are public.
Sim.DiffProc documentation built on Nov. 8, 2020, 4:27 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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 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:
dataa univariate time series (tsobject).driftanddiffusionindicate drift and diffusion coefficient of the SDE, is anexpressionof two variablest,xandthetanames of the parameters, and must be nominated by a vector oftheta = (theta[1], theta[2],..., theta[p])for reasons of symbolic derived in approximation methods.startmust 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 thedriftanddiffusioncoefficient.- The
pmleargument must be acharacterstring 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.methodselect the optimization method ("L-BFGS-B"is used by default), and further arguments to pass tooptimfunction.lowerandupperbounds on the variables for theBrentorL-BFGS-Bmethod.
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 byfitsde.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:
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".
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:
coef(fitmod)
We can use the summary function to produce result summaries of output object:
summary(fitmod)
vcov for variance-covariance matrice, and extract log-likelihood by logLik:
vcov(fitmod) logLik(fitmod) AIC(fitmod) BIC(fitmod)
Computes confidence intervals for one or more parameters in a fitted SDE:
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:
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:
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:
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:
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":
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:
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":
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:
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:
## 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:
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, 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.
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:
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:
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
snssdekd()&dsdekd()&rsdekd()- Monte-Carlo Simulation and Analysis of Stochastic Differential Equations.bridgesdekd()&dsdekd()&rsdekd()- Constructs and Analysis of Bridges Stochastic Differential Equations.fptsdekd()&dfptsdekd()- Monte-Carlo Simulation and Kernel Density Estimation of First passage time.MCM.sde()&MEM.sde()- Parallel Monte-Carlo and Moment Equations for SDEs.TEX.sde()- Converting Sim.DiffProc Objects to LaTeX.fitsde()- Parametric Estimation of 1-D Stochastic Differential Equation.
References
-
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.
-
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
-
Iacus SM (2008). Simulation and Inference for Stochastic Differential Equations: With R Examples. Springer-Verlag, New York.
Try the Sim.DiffProc package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.