[Yuima-commits] r200 - pkg/yuimadocs/inst/doc/JSS
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jul 25 11:42:03 CEST 2012
Author: iacus
Date: 2012-07-25 11:42:03 +0200 (Wed, 25 Jul 2012)
New Revision: 200
Added:
pkg/yuimadocs/inst/doc/JSS/article-new.Rnw
Modified:
pkg/yuimadocs/inst/doc/JSS/article.Rnw
Log:
revision
Added: pkg/yuimadocs/inst/doc/JSS/article-new.Rnw
===================================================================
--- pkg/yuimadocs/inst/doc/JSS/article-new.Rnw (rev 0)
+++ pkg/yuimadocs/inst/doc/JSS/article-new.Rnw 2012-07-25 09:42:03 UTC (rev 200)
@@ -0,0 +1,1212 @@
+\documentclass[article]{jss}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% declarations for jss.cls %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%\usepackage{Sweave}
+%\usepackage{Rd}
+
+\usepackage{amsmath,amssymb}
+
+\usepackage[utf8x]{inputenc}
+
+%% almost as usual
+\author{Alexandre Brouste\\University of Le Mans \And
+ Masaaki Fukasawa\\Osaka University\And
+ Hideitsu Hino\\Waseda University\And
+ Stefano M. Iacus\\University of Milan \AND
+ Kengo Kamatani\\The University of Tokyo \And
+ Yuta Koike\\The University of Tokyo\And
+ Hiroki Masuda\\Kyushu University \And
+ Ryosuke Nomura\\The University of Tokyo\AND
+ Yasutaka Shimuzu\\Osaka University \And
+ Masayuki Uchida\\Osaka University \And
+ Nakahiro Yoshida\\The University of Tokyo
+ }
+\title{The YUIMA Project: a Computational Framework for Simulation and Inference of Stochastic Differential Equations}
+
+%% for pretty printing and a nice hypersummary also set:
+\Plainauthor{Alexandre Brouste, Stefano M. Iacus, Hiroki Masuda, Masayuki Uchida, Nakahiro Yoshida} %% comma-separated
+\Plaintitle{YUIMA: Simulation and Inference for SDE} %% without formatting
+\Shorttitle{The YUIMA Project} %% a short title (if necessary)
+
+%% an abstract and keywords
+\Abstract{
+The Yuima Project is an open source and collaborative effort aimed at developing the \proglang{R} package named \pkg{yuima} for simulation and inference of stochastic differential equations.
+
+In the \pkg{yuima} package stochastic differential equations can be of very abstract type, multidimensional, driven by Wiener process or fractional Brownian motion with general Hurst parameter, with or without jumps specified as L\'evy noise.
+
+The \pkg{yuima} package is intended to offer the basic infrastructure on which complex models and inference procedures can be built on. This paper explains the design of the \pkg{yuima} package and provides some examples of applications.
+}
+\Keywords{inference for stochastic processes, simulation, stochastic differential equations}
+\Plainkeywords{inference for stochastic processes, simulation, stochastic differential equations} %% without formatting
+%% at least one keyword must be supplied
+
+%% publication information
+%% NOTE: Typically, this can be left commented and will be filled out by the technical editor
+%% \Volume{13}
+%% \Issue{9}
+%% \Month{September}
+%% \Year{2004}
+%% \Submitdate{2004-09-29}
+%% \Acceptdate{2004-09-29}
+
+%% The address of (at least) one author should be given
+%% in the following format:
+\Address{
+ Stefano M. Iacus\\
+ Department of Economics, Management and Quantitative Methods\\
+ University of Milan\\
+ Via Conservatorio 7, 20122 Milan, Italy\\
+ E-mail: \email{stefano.iacus at unimi.it}\\
+ URL: \url{http://www.economia.unimi.it/~iacus/}
+\\
+\\
+ Nakahiro Yoshida\\
+ Graduate School of Mathematical Science\\
+ University of Tokyo\\
+ 3-8-1 Komaba, Meguro-ku, Tokyo 153-8914, Japan\\
+ E-mail: \email{nakahiro at ms.u-tokyo.ac.jp}\\
+ URL: \url{http://www.ms.u-tokyo.ac.jp/~nakahiro/hp-naka-e}
+}
+%% It is also possible to add a telephone and fax number
+%% before the e-mail in the following format:
+%% Telephone: +43/1/31336-5053
+%% Fax: +43/1/31336-734
+
+
+
+\SweaveOpts{prefix.string=yuima}
+\SweaveOpts{echo=TRUE}
+
+%% preamble for Rnw files
+<<print=FALSE, echo=FALSE>>=
+options(prompt="R> ")
+options(width=60)
+@
+
+\def\ve{{\varepsilon}}
+\def\de{{\rm d}}
+\def\dE{{\mathbb E}}
+\def\diag{\mathop{\rm diag}\nolimits}
+
+%% for those who use Sweave please include the following line (with % symbols):
+%% need no \usepackage{Sweave.sty}
+%% end of declarations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+\begin{document}
+
+%% include your article here, just as usual
+%% Note that you should use the \pkg{}, \proglang{} and \code{} commands.
+
+%\section[About Java]{About \proglang{Java}}
+%% Note: If there is markup in \(sub)section, then it has to be escape as above.
+
+
+
+
+
+
+
+
+\section{Introduction}
+The plan of the YUIMA Project is to define the bases for a complete environment for simulation and inference for stochastic processes via an \proglang{R} \citep{ERRE} package called \pkg{yuima}.
+ The package \pkg{yuima} provides
+an object-oriented programming environment
+for simulation and statistical inference
+for stochastic processes by \proglang{R}.
+The \pkg{yuima} package adopts the S4 system of classes and methods \citep{chambers98}.
+
+Under this framework,
+the \pkg{yuima} package also supplies various functions
+to execute simulation and statistical analysis.
+Both categories of procedures may depend \textcolor{red}{on} each other.
+Statistical inference often requires a simulation technique
+as a subroutine, and a certain simulation method
+needs to fix a tuning parameter by applying
+a statistical methodology.
+It is especially the case of stochastic processes
+because most of expected values involved
+do not admit an explicit expression.
+The \pkg{yuima} package facilitates comprehensive, systematic approaches
+to the solution.
+
+
+Stochastic differential equations are
+commonly used
+to model random evolution along continuous or
+practically continuous time, such as
+the random movements of a stock price.
+Theory of statistical inference for
+stochastic differential equations already has a fairly long history,
+more than three decades, but it is still developing quickly new
+methodologies and expanding the area.
+The formulas produced by the theory are usually very sophisticated,
+which makes it difficult for standard users not necessarily
+familiar with this field to enjoy utilities.
+\textcolor{red}{For example, by taking advantage of the analytic approach, the asymptotic expansion method for computing
+option prices (i.e., expectation of an irregular functional of
+a stochastic process) provides precise approximation values
+instantaneously. The expansion formula, which has a long expression involving more than 900 terms of multiple integrals, is already coded in the \pkg{yuima} package for generic diffusion processes.}
+
+
+The \pkg{yuima} package delivers up-to-date methods as a package
+onto the desk of the user working
+with simulation and/or statistics for stochastic differential equations.
+In the \pkg{yuima} package stochastic differential equations can be of very abstract type,
+multidimensional, driven by Wiener process or fractional Brownian motion
+with general Hurst parameter, with or without jumps specified as L\'evy noise.
+
+
+
+The \pkg{yuima} package is intended to offer the basic infrastructure on which complex models and inference procedures can be built on. This paper explains the design of the \pkg{yuima} package and provides some examples of applications.
+The paper is organised as follows. Section \ref{sec2} is an overview about the package. Section \ref{sec3} describe the way models are specified in \pkg{yuima}. Section \ref{sec4} explains asymptotic expansion methods. Section \ref{sec5} is a review of basic inference procedures. Finally, Section \ref{sec6} explains additional details and the roadmap of the YUIMA Project.
+
+
+Although we assume that the reader of this paper has a basic knowledge of the \proglang{R} language, most of the examples are easy to be understood by anyone.
+
+\section{The \pkg{yuima} Package}\label{sec2}
+The package \pkg{yuima} depends on some other packages, like \pkg{zoo} \citep{zoo}, which can be installed separately.
+The package \pkg{zoo} is used internally to store time series data. This dependence may change in the future adopting a more flexible class for internal storage of time series.
+
+\subsection{How to Obtain the Package}
+The \pkg{yuima} package is hosted on R-Forge and the web page
+of the Project is \url{http://r-forge.r-project.org/projects/yuima}.
+The R-Forge page contains the latest development version, and stable
+version of the package will also be available through CRAN.
+Development versions of the package are not supposed to be stable or functional, thus
+the occasional user should consider to install the stable version first.
+The package can be installed from R-Forge using
+\code{install.packages("yuima",repos="http://R-Forge.R-project.org")}.
+
+
+\subsection{The Main Object and Classes}
+Before discussing the methods for simulation and inference for stochastic processes solutions to stochastic differential equations, here we discuss the main classes in the package.
+As mentioned there are different classes of objects defined in the \pkg{yuima} package and
+the main class is called the \code{yuima-class}. This class is composed of several slots.
+Figure \ref{fig:classes} represents the different classes and their slots.
+\begin{figure}[th]
+\centering{\includegraphics[width=\textwidth]{yuima-class}}
+\caption{The main classes in the \pkg{yuima} package.}
+\label{fig:classes}
+\end{figure}
+The different slots do not need to be all present at the same time. For example, in case one wants to simulate a stochastic process, only the slots \code{model} and \code{sampling} should be present, while the slot \code{data} will be filled by the simulator.
+We now discuss in details the \textcolor{red}{different objects} separately.
+
+
+
+\subsection{The \code{yuima.model} Class}
+At present, in \pkg{yuima} three main classes of stochastic differential equations can be easily specified. All multidimensional and eventually as parametric models.
+\begin{itemize}
+\item diffusions $\displaystyle \de X_t=a(t,X_t)dt + b(t,X_t)\de W_t $, where $W_t$ is a standard Brownian motion;
+\item fractional Gaussian noise, with $H$ the Hurst parameter
+$$ \de X_t=a(t,X_t)dt + b(t,X_t)\de W_t^{H}; $$
+\item diffusions with jumps and L\'evy processes solution to
+$$
+\begin{aligned}
+\de X_t = & \,\,\, a(X_t)\de t + b(X_t)\de W_t + \int\limits_{|z|>1}\!\!\! c(X_{t-},z)\mu(\de t,\de z) \\
+&{}+\!\! \int\limits_{0<|z|\le 1}\!\!\! c(X_{t-},z)\{\mu(\de t,\de z) - \nu(\de z)\de t\}.
+\end{aligned}
+$$
+\end{itemize}
+
+The \code{yuima.model} class contains information about the stochastic differential equation of interest. The constructor function \code{setModel} is
+used to give a mathematical description of the stochastic differential equation.
+All functions in the package are assumed to get as much information
+as possible from the model instead of replicating the same code everywhere.
+If there are missing pieces of information, we may change or extend the description
+of the model.
+\\
+\\
+An object of class \code{yuima.model} contains several slots listed below.
+To see inside its structure, one can use the \proglang{R} command \code{str} on a \pkg{yuima} object.
+\begin{itemize}
+\item \code{drift} is an \proglang{R} vector of expressions which contains the drift specification;
+\item \code{diffusion} is itself a list of 1 slot which describes the diffusion
+ coefficient matrix relative to first noise;
+\item \code{hurst} is the Hurst index of the fractional Brownian motion, by default \code{0.5} meaning a standard Brownian motion;
+\item \code{parameter} which is a short name for ``parameters'' which
+ is a list with the following entries:
+\begin{itemize}
+\item \code{all} contains the names of all the parameters found
+ in the diffusion and drift coefficient;
+\item \code{common} contains the names of the parameters in common between the drift and diffusion coefficients;
+\item \code{diffusion} contains the parameters belonging to the diffusion coefficient;
+\item \code{drift} contains the parameters belonging to the drift coefficient;
+\item \code{jump} contains the parameters belonging to the coefficient of the L\'evy noise;
+\item \code{measure} contains the parameters \textcolor{red}{describing the L\'evy} measure;
+\end{itemize}
+\item \code{measure} specifies the measure of the L\'evy noise;
+\item \code{measure.type} a switch to specify how the L\'evy measure is described;
+\item \code{state.variable} and \code{time.variable}, by default,
+ are assumed to be \code{x} and \code{t} but the user can freely choose them.
+ The \code{yuima.model} class assumes that the user either \textcolor{red}{uses} default
+ names for \code{state.variable} and \code{time.variable} variables or specify his own names. All
+ the rest of the symbols are considered parameters and distributed accordingly
+ in the \code{parameter} slot.
+\item \code{jump.variable} the name of the variable used in the description of the L\'evy component;
+\item \code{solve.variable} contains a vector of variable names, each element corresponds to the
+ name of the solution variable (left-hand-side) of each equation in the model, in the corresponding order.
+\item \code{noise.number} indicates the number of sources of noise.
+\item \code{xinit} initial value of the stochastic differential equation;
+\item \code{equation.number} represents the number of equations, i.e., the number of one dimensional
+ stochastic differential equations.
+\item \code{dimension} reports the dimensions of the parameter space. It is a list of the
+ same length of \code{parameter} with the same names.
+\item \code{J.flag} for internal use only.
+\end{itemize}
+As seen in the above, the parameter space is accurately described internally in a \code{yuima} object because in inference for stochastic differential equations, estimators of different parameters have different properties. Usually, the rate of convergence for estimators in the diffusion coefficient are similar to the ones in the i.i.d. sampling while estimators of parameters in the drift coefficient are slower or, in some cases, not even consistent. The \pkg{yuima} always tries to implement the best statistical inference for the given model under the observed sampling scheme.
+
+\section{Model Specification}\label{sec3}
+In order to show \textcolor{red}{how general the approach is} in the \pkg{yuima} package we present some examples.
+
+\subsection{Diffusion Processes}
+Assume that we want to describe the following stochastic differential equation
+$$\de X_t = -3 X_t \de t + \frac{1}{1+X_t^2}\de W_t$$
+This is done in \pkg{yuima} \textcolor{red}{by} specifying the drift and diffusion coefficients as plain mathematical expressions
+<<print=FALSE,echo=FALSE,results=hide>>=
+library(yuima)
+@
+<<echo=TRUE, print=FALSE,results=hide>>=
+mod1 <- setModel(drift = "-3*x", diffusion = "1/(1+x^2)")
+@
+At this point, the package fills the proper slots of the \code{yuima} object
+<<>>=
+str(mod1)
+@
+
+For the above, it is possible to see that the jump coefficient is void and the Hurst parameter is set to 0.5, because this corresponds to the standard Brownian motion.
+Now, with \code{mod1} in hands, it is very easy to simulate a trajectory of the process as follows
+<<echo=TRUE, print=FALSE,fig=TRUE,width=9,height=4,results=hide>>=
+set.seed(123)
+X <- simulate(mod1)
+plot(X)
+@
+
+\noindent
+The \code{simulate} function fills in addition the two slots \code{data} and \code{sampling} of the \code{yuima} object.
+<<>>=
+str(X,vec.len=2)
+@
+
+\subsection{Specification of Parametric Models}
+When a parametric model like
+$$\de X_t = -\theta X_t \de t + \frac{1}{1+X_t^\gamma}\de W_t$$
+is specified, \pkg{yuima} attempts to distinguish the parameters' names from the ones of the state and time variables
+<<echo=TRUE, print=FALSE,results=hide>>=
+mod2 <- setModel(drift = "-theta*x", diffusion = "1/(1+x^gamma)")
+@
+<<>>=
+str(mod2)
+@
+In order to simulate the parametric model it is necessary to specify the values of the parameters as the next code shows
+<<echo=TRUE, print=FALSE,fig=TRUE,height=4,results=hide>>=
+set.seed(123)
+X <- simulate(mod2,true.param=list(theta=1,gamma=3))
+plot(X)
+@
+
+\subsection{Multidimensional Processes}
+Next is an example with two stochastic differential equations driven by three independent Brownian motions
+$$
+\begin{aligned}
+\de X^1_t &= -3 X^1_t \de t + \de W^1_t + X^2_t \de W^3_t\\
+\de X^2_t &= -(X^1_t + 2 X^2_t) \de t + X^1_t \de W^1_t + 3 \de W^2_t
+\end{aligned}
+$$
+but this has to be organized into matrix form
+$$
+\left(\begin{array}{c}\de X^1_t \\\de X^2_t\end{array}\right)=
+\left(\begin{array}{c} -3 X^1_t \\ -X^1_t - 2X^2_t\end{array}\right)\de t +
+\left[\begin{array}{ccc}1 & 0 & X^2_t \\X^1_t & 3 & 0\end{array}\right]
+\left(\begin{array}{c}\de W_t^1 \\ \de W_t^2 \\ \de W_t^3\end{array}\right)
+$$
+<<echo=TRUE, print=FALSE,results=hide>>=
+sol <- c("x1","x2") # variable for numerical solution
+a <- c("-3*x1","-x1-2*x2") # drift vector
+b <- matrix(c("1","x1","0","3","x2","0"),2,3) # diffusion matrix
+mod3 <- setModel(drift = a, diffusion = b, solve.variable = sol)
+@
+Again, this model can be easily simulated
+<<echo=TRUE, print=FALSE,fig=TRUE,width=9,height=4,results=hide>>=
+set.seed(123)
+X <- simulate(mod3)
+plot(X, plot.type="single",lty=1:2)
+@
+
+\noindent
+But it is also possible to specify more complex models like the following
+$$
+ \begin{cases}
+ &\de X_t^1 = X_t^2 \left| X_t^1 \right|^{2/3} \de W_t^1, \\
+ &\de X_t^2 = g(t)\de X_t^3, \\
+ &\de X_t^3 = X_t^3( \mu \de t + \sigma (\rho \de W_t^1 + \sqrt{1-\rho^2} \de W_t^2))
+ \end{cases},
+$$
+where $g(t) = 0.4 + (0.1 + 0.2t) e^{-2t}$.
+
+\subsection{Fractional Gaussian Noise}
+In order to specify a stochastic differential equation driven by fractional Gaussian noise it is necessary to specify the value of the Hurst parameter. For example, if we want to specify the following model
+$$\de Y_t = 3 Y_t \de t + \de W_t^H$$
+we proceed as follows
+<<echo=TRUE, print=FALSE,fig=TRUE,width=9,height=4,results=hide>>=
+mod4 <- setModel(drift="3*y", diffusion=1, hurst=0.3, solve.var="y")
+set.seed(123)
+X <- simulate(mod4,sampling=setSampling(n=1000))
+plot(X)
+@
+
+\noindent
+In this case, the appropriate slot is now filled
+@
+<<>>=
+str(mod4)
+@
+The user can choose \textcolor{red}{between two} simulation schemes, namely the Cholesky method and \citet{WoodChan} method.
+
+\subsection{L\'evy Processes}
+Jump processes can be specified in different ways in mathematics and hence in \pkg{yuima} package.
+Let $Z_t$ be a Compound Poisson Process (i.e., jumps size follow some distribution, like the Gaussian law and jumps occur at Poisson times).
+Then it is possible to consider the following SDE which involves jumps
+$$\de X_t = a(X_t)\de t + b(X_t)\de W_t + \de Z_t$$
+In the next example we consider a compound Poisson process with intensity $\lambda=10$ with Gaussian jumps.
+This model can be specified in \code{setModel} using the argument \code{measure.type="CP"}
+A simple Ornstein-Uhlenbeck process with Gaussian jumps
+$$\de X_t = -\theta X_t \de t + \sigma \de W_t + \de Z_t$$
+is specified as
+<<echo=TRUE, print=FALSE,fig=TRUE,width=9,height=4,results=hide>>=
+mod5 <- setModel(drift=c("-theta*x"), diffusion="sigma",
+ jump.coeff="1", measure=list(intensity="10", df=list("dnorm(z, 0, 1)")),
+ measure.type="CP", solve.variable="x")
+set.seed(123)
+X <- simulate(mod5, true.p=list(theta=1,sigma=3),sampling=setSampling(n=1000))
+plot(X)
+@
+
+\noindent
+Another possibility is to specify the jump component via the L\'evy measure. Without going into too much details, here is an example of specification of a simple Ornstein-Uhlenbeck process with IG (Inverse Gaussian) L\'evy measure
+$$\de X_t = -x \de t + \de Z_t$$
+<<echo=TRUE, print=FALSE,fig=TRUE,width=9,height=4,results=hide>>=
+mod6 <- setModel(drift="-x", xinit=1, jump.coeff="1",
+ measure.type="code", measure=list(df="rIG(z, 1, 0.1)"))
+set.seed(123)
+X <- simulate(mod6, sampling=setSampling(Terminal=10, n=10000))
+plot(X)
+@
+
+\subsection{Specification of Generic Models}
+In general, the \pkg{yuima} package allows to specify a large family of models solutions to
+$$\de X_t =a(X_t)\de t + b(X_t)\de W_t + c(X_t)\de Z_t$$
+using the following interface
+<<echo=TRUE, print=FALSE,eval=FALSE>>=
+setModel(drift, diffusion, hurst = 0.5, jump.coeff,
+ measure, measure.type, state.variable = "x",
+ jump.variable = "z", time.variable = "t", solve.variable, xinit)
+@
+The \pkg{yuima} package implements many multivariate Random Numbers Generators (RNG) which are needed to simulate L\'evy paths including
+\code{rIG} (Inverse Gaussian), \code{rNIG} (Normal Inverse Gaussian), \code{rbgamma} (Bilateral Gamma), \code{rngamma} (Normal Gamma) and \texttt{rstable} (Stable Laws).
+Other user-defined RNG can be used freely.
+
+\subsection{Simulation, Sampling and Subsampling}
+The \code{simulate} function simulates \code{yuima} models according to Euler-Maruyama scheme in the presence of non-fractional diffusion noise and L\'evy jumps and uses the Cholesky or the Wood and Chan methods for the fractional Gaussian noise.
+The \code{simulate} function accepts several arguments including the description of the sampling structure, which is an object of type \code{yuima.sampling}. The \code{setSampling} \textcolor{red}{allows} for the specification of different sampling parameters including random sampling. Further, the \code{subsampling} \textcolor{red}{allows us} to subsample a trajectory of a simulated stochastic differential equation or a given time series in the \code{yuima.data} slot of a \code{yuima} object.
+Sampling and subsampling can be specified jointly as arguments to the \code{simulate} function. This is convenient if one wants to simulate data at very high frequency but then return only low frequency data for inference or other applications. \textcolor{red}{In what follows we explain how to specify arguments of these \pkg{yuima} functions. Complete details can be found in the man pages of the \pkg{yuima} package.}
+
+Assume that we want to simulate this model
+$$
+\begin{aligned}
+\de X^1_t &= -\theta X^1_t \de t + \de W^1_t + X^2_t \de W^3_t\\
+\de X^2_t &= -(X^1_t + \gamma X^2_t) \de t + X^1_t \de W^1_t + \delta \de W^2_t
+\end{aligned}
+$$
+Now we prepare the model using the \code{setModel} constructor function
+<<setModel1>>=
+sol <- c("x1","x2") # variable for numerical solution
+b <- c("-theta*x1","-x1-gamma*x2") # drift vector
+s <- matrix(c("1","x1","0","delta","x2","0"),2,3) # diffusion matrix
+mymod <- setModel(drift = b, diffusion = s, solve.variable = sol)
+@
+Suppose now that we want to simulate the process on a regular grid on the interval $[0,3]$ and $n=3000$ observations. We can prepare the sampling structure as follows
+<<>>=
+samp <- setSampling(Terminal=3, n=3000)
+@
+and look inside it
+<<>>=
+str(samp)
+@
+As seen from the output, the sampling structure is \textcolor{red}{quite rich} and we will show how to specify \textcolor{red}{some} of the slots \textcolor{red}{later on}.
+We now simulate this process specifying the \code{sampling} argument to \code{simulate}
+<<>>=
+set.seed(123)
+X2 <- simulate(mymod, sampling=samp)
+@
+now the sampling structure is recorded along with the data in the \code{yuima} object \code{X2}
+<<>>=
+str(X2 at sampling)
+@
+\subsubsection{Subsampling}
+The sampling structure can be used to operate subsampling. Next example shows how to perform Poisson random sampling, with two independent Poisson processes, one per coordinate of \code{X2}.
+<<sub1>>=
+newsamp <- setSampling(
+random=list(rdist=c( function(x) rexp(x, rate=10),
+function(x) rexp(x, rate=20))) )
+str(newsamp)
+@
+In the above we have specified two independent exponential distributions to represent \textcolor{red}{Poisson arrival times}. Notice that the slot \code{regular} is now set to \code{FALSE}.
+Now we subsample the original trajectory of \code{X2} using the \code{subsampling} function
+<<sub2,fig=TRUE,width=9,height=5>>=
+newdata <- subsampling(X2, sampling=newsamp)
+plot(X2,plot.type="single", lty=c(1,3),ylab="X2")
+points(get.zoo.data(newdata)[[1]],col="red")
+points(get.zoo.data(newdata)[[2]],col="green",pch=18)
+@
+
+Or we can operate a deterministic sampling specifying two different regular frequencies
+<<sub3,fig=TRUE,width=9,height=5>>=
+newsamp <- setSampling(delta=c(0.1,0.2))
+newdata <- subsampling(X2, sampling=newsamp)
+plot(X2,plot.type="single", lty=c(1,3),ylab="X2")
+points(get.zoo.data(newdata)[[1]],col="red")
+points(get.zoo.data(newdata)[[2]],col="green", pch=18)
+@
+
+Again one can look at the structure of the sampling structure.
+
+Subsampling can be used within the \code{simulate} function. What is usually done in simulation studies, is to simulate the process at very high frequency but then use data for estimation at a lower frequency.
+This can be done in a single step in the following way.
+<<sub4,fig=TRUE,width=9,height=5>>=
+set.seed(123)
+Y.sub <- simulate(mymod,sampling=setSampling(delta=0.001,n=1000), subsampling=setSampling(delta=0.01,n=100))
+set.seed(123)
+Y <- simulate(mymod, sampling=setSampling(delta=0.001,n=1000))
+plot(Y, plot.type="single")
+points(get.zoo.data(Y.sub)[[1]],col="red")
+points(get.zoo.data(Y.sub)[[2]],col="green",pch=18)
+@
+
+In the previous code we have simulated the process twice just to show the effect of the subsampling but the reader should use only the line which outputs the simulation to \code{Y.sub} as seen in the next plot.
+<<sub5,fig=TRUE,width=9,height=5>>=
+plot(Y.sub, plot.type="single")
+@
+
+
+
+\section{Asymptotic Expansion}\label{sec4}
+The \pkg{yuima} package can handle asymptotic expansion of functionals of $d$-dimensional diffusion process
+$$\de X_t^\ve = a(X_t^\ve,\ve)\de t + b(X_t^\ve,\ve)\de W_t, \qquad \ve \in(0,1]$$
+with $W_t$ \textcolor{red}{an} $r$-dimensional Wiener process, i.e., $W_t=(W_t^1, \ldots, W_t^r)$.
+The functional is expressed in the following abstract form
+\begin{equation}
+\label{yuima:func1}
+F^\ve(X_t^\ve) = \sum_{\alpha=0}^r \int_0^T f_\alpha(X_t^\ve,\de)\de W_t^\alpha + F(X_t^\ve,\ve), \qquad W_t^0=t.
+\end{equation}
+A typical example of application is the case of Asian option pricing. For example, in the Black \& Scholes model
+\begin{equation}
+\label{eqyuima:gbm1}
+\de X_t^\ve = \mu X_t^\ve\de t + \ve X_t^\ve \de W_t
+\end{equation}
+the price of the option is of the form
+$$ \dE\left\{ \max\left( \frac{1}{T}\int_0^T X_t^\ve \de t - K,0\right)\right\}.$$
+Thus the functional of interest is
+$$F^\ve(X_t^\ve) = \frac{1}{T}\int_0^T X_t^\ve \de t, \qquad r=1$$
+with
+$$
+f_0(x,\ve) = \frac{x}{T}, \quad f_1(x,\ve)=0, \quad F(x, \ve) = 0
+$$
+in \eqref{yuima:func1}.
+So, the call option price requires the composition of a smooth functional
+$$F^\ve(X_t^\ve) = \frac{1}{T}\int_0^T X_t^\ve \de t, \qquad r=1$$
+with the irregular function
+$$
+\max(x-K,0)
+$$
+Monte Carlo methods require a huge number of simulations to get the desired accuracy of the
+calculation of the price, while asymptotic expansion of $F^\ve$ provides very accurate approximations.
+The \pkg{yuima} package provides functions to construct the functional $F^\ve$, and automatic asymptotic expansion based on Malliavin calculus \citep{Yoshida92b} starting from a \pkg{yuima} object.
+Here is an example. Consider a simple geometric Brownian motion of equation \eqref{eqyuima:gbm1} with $\mu=1$ and $X_0=1$. We define the model and the functional below:
+<<echo=TRUE, print=FALSE, results=hide>>=
+model <- setModel(drift = "x", diffusion = matrix( "x*e", 1,1))
+T <- 1
+xinit <- 1
+K <- 1
+f <- list( expression(x/T), expression(0))
+F <- 0
+e <- 0.5
+yuima <- setYuima(model = model, sampling = setSampling(Terminal=T, n=1000))
+yuima <- setFunctional( yuima, f=f,F=F, xinit=xinit,e=e)
+@
+this time the \code{setFunctional} command fills the appropriate slots inside the \code{yuima} object
+<<echo=TRUE, print=FALSE>>=
+str(yuima at functional)
+@
+Then, to obtain the first term in the asymptotic expansion (i.e., for $\ve=0$), it is as easy as calling the function \code{F0} on the \code{yuima} object:
+<<echo=TRUE, print=FALSE>>=
+F0 <- F0(yuima)
+F0
+@
+so the option price approximation is
+<<>>=
+max(F0-K,0) # asian call option price
+@
+We can go up to the first order approximation adding one term to the expansion
+<<echo=TRUE, print=FALSE,results=hide>>=
+rho <- expression(0)
+epsilon <- e # noise level
+g <- function(x) {
+ tmp <- (F0 - K) + (epsilon * x)
+ tmp[(epsilon * x) < (K-F0)] <- 0
+ tmp
+}
+asymp <- asymptotic_term(yuima, block=10, rho, g)
+@
+and the final value is
+<<echo=TRUE, print=FALSE>>=
+asymp$d0 + e * asymp$d1 # asymp. exp. of asian call price
+@
+
+
+
+\section{Inference for Stochastic Processes}\label{sec5}
+The \pkg{yuima} implements several optimal techniques for parametric, semi- and non-parametric estimation of (multidimensional) stochastic differential equations.
+Although most of the examples in this section are given on simulated data, the main way to fill up the \code{data} slot of a \code{yuima} object is to use the function \code{setYuima}. The function \code{setYuima} sets various slots of the \code{yuima} object. In particular, to estimate a \code{yuima.model} called \code{mod} on the data \code{X} one can use a code like this \code{my.yuima <- setYuima(data=setData(X), model=mod)} and then pass \code{my.yuima} to the inference functions as described in what follows.
+
+\subsection{Quasi Maximum Likelihood Estimation}
+Consider the multidimensional diffusion process
+$$
+\de X_t = b(\theta_2,X_t)\de t + \sigma(\theta_1, X_t) \de W_t
+$$
+where $W_t$ is an {$r$}-dimensional standard Wiener process
+independent of the initial value $X_0=x_0$.
+Quasi-MLE assumes the following approximation of the true log-likelihood for multidimensional diffusions
+{\small
+\begin{eqnarray}\label{qlik}
+\ell_n({\bf X}_n,\theta)
+=-\frac12\sum_{i=1}^n\left\{\log\text{det}(\Sigma_{i-1}(\theta_1))
++\frac{1}{\Delta_n}
+\Sigma_{i-1}^{-1}(\theta_1)[\Delta X_i-\Delta_n b_{i-1}(\theta_2)]^{\otimes 2}\right\}
+\end{eqnarray}}
+where $\theta=(\theta_1, \theta_2)$, $\Delta X_i=X_{t_i}-X_{t_{i-1}}$, $\Sigma_i(\theta_1)=\Sigma(\theta_1,X_{t_i})$, $b_i(\theta_2)=b(\theta_2,X_{t_i})$, $\Sigma=\sigma^{\otimes 2}$, \textcolor{red}{$A^{\otimes 2}= A A^T$} and $A^{-1}$ the inverse of $A$, $A[B]^{\otimes 2} = B^T A B$. Then, \citep[see e.g.,][]{Yoshida92, Kessler97}, the QML estimator of $\theta$ is
+\textcolor{red}{$$\hat\theta_n=\arg\max_\theta \ell_n({\bf X}_n,\theta)$$}
+
+\textcolor{red}{The \pkg{yuima} package implements QML estimation via the \code{qmle} function. The interface and the output of the \code{qmle} function are made as similar as possible to those of the standard \code{mle} function in the \pkg{stats4} package of the basic \proglang{R} system. The main arguments to \code{qmle} consist of a \code{yuima} object and initial values (\code{start}) for the optimizer. The \code{yuima} object must contain the slots \code{model} and \code{data}. The \code{start} argument 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 \code{yuima} object. Optionally, one can specify named lists of \code{upper} and \code{lower} bounds
+to identify the search region of the optimizer. The standard optimizer is \code{BFGS} when no bounds are specified. If bounds are specified then \code{L-BFGS-B} is used. More optimizers can be added in the future.}
+
+
+
+As an example, we consider the simple model
+\begin{equation}
+\de X_t = (2-\theta_2 X_t) \de t + (1+X_t^2)^{\theta_1} \de W_t
+\label{eq:model1}
+\end{equation}
+with $\theta_1=0.2$ and $\theta_2=0.3$. We generate sampled data $X_{t_i}$, with $t_i = i \cdot n^{-\frac23}$.
+<<echo=TRUE, print=FALSE,results=hide>>=
+ymodel <- setModel(drift="(2-theta2*x)", diffusion="(1+x^2)^theta1")
+n <- 750
+ysamp <- setSampling(Terminal = n^(1/3), n = n)
+yuima <- setYuima(model = ymodel, sampling = ysamp)
+set.seed(123)
+yuima <- simulate(yuima, xinit = 1, true.parameter = list(theta1 = 0.2, theta2 = 0.3))
+@
+\textcolor{red}{Now the \code{yuima} object contains both the \code{model} and \code{data} slots. We set the initial values for the optimizer as $\theta_1=\theta_2=0.5$ and we specify them as a named list called \code{param.init}. We can now use the function \code{qmle} to estimate the parameters $\theta_1$ and $\theta_2$ as follows}
+<<echo=TRUE, print=FALSE,results=hide>>=
+param.init <- list(theta2=0.5,theta1=0.5)
+upp.par <- list(theta1=0, theta2=0)
+low.par <- list(theta1=1, theta2=1)
+mle1 <- qmle(yuima, start = param.init, lower = upp.par, upper = low.par)
+@
+\textcolor{red}{where \code{upp.par} and \code{low.par} are the upper and lower bounds of the search region used by the optimizer (\code{L-BFGS-B} in this case).}
+\textcolor{red}{The} estimated coefficients are extracted from the output object \code{mle1} as follows
+<<echo=TRUE, print=FALSE>>=
+summary(mle1)
+@
+%Notice the interface and the output of the \code{qmle} is quite similar to the ones of the standard \code{mle} function of the \pkg{stats4} package of the base \proglang{R} system.
+
+\subsection{Adaptive Bayes Estimation}
+Consider again the diffusion process solution to
+\begin{equation}
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 200
More information about the Yuima-commits
mailing list