[Yuima-commits] r11 - in pkg: . yuimadocs yuimadocs/inst yuimadocs/inst/doc
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Nov 7 19:50:57 CET 2009
Author: iacus
Date: 2009-11-07 19:50:57 +0100 (Sat, 07 Nov 2009)
New Revision: 11
Added:
pkg/yuimadocs/
pkg/yuimadocs/DESCRIPTION
pkg/yuimadocs/inst/
pkg/yuimadocs/inst/doc/
pkg/yuimadocs/inst/doc/yuima.Rnw
pkg/yuimadocs/inst/doc/yuima_1.Rnw
pkg/yuimadocs/inst/doc/yuima_2.Rnw
pkg/yuimadocs/inst/doc/yuima_3.Rnw
pkg/yuimadocs/inst/doc/yuima_4.Rnw
Log:
package skeleton for yuimadoc package, not yet working
Added: pkg/yuimadocs/DESCRIPTION
===================================================================
--- pkg/yuimadocs/DESCRIPTION (rev 0)
+++ pkg/yuimadocs/DESCRIPTION 2009-11-07 18:50:57 UTC (rev 11)
@@ -0,0 +1,11 @@
+Package: yuimadocs
+Type: Package
+Title: The YUIMA Project package Doc Collection
+Version: 0.0.1
+Date: 2009-11-07
+Depends: yuima, methods, zoo
+Author: Yuima project team.
+Maintainer: Who to complain to <yourfault at somewhere.net>
+Description: package for the YUIMA Project. The YUIMA Project Team is: Yoshida, N., Iacus, S.M., Hino, H.,
+License: GPL (>= 2)
+
Property changes on: pkg/yuimadocs/DESCRIPTION
___________________________________________________________________
Name: svn:executable
+ *
Added: pkg/yuimadocs/inst/doc/yuima.Rnw
===================================================================
--- pkg/yuimadocs/inst/doc/yuima.Rnw (rev 0)
+++ pkg/yuimadocs/inst/doc/yuima.Rnw 2009-11-07 18:50:57 UTC (rev 11)
@@ -0,0 +1,46 @@
+\documentclass[a4paper]{article}
+
+\title{yuimaS4 package documentation v1.0}
+\author{}
+
+\SweaveOpts{echo=FALSE}
+\usepackage{a4wide}
+\usepackage{Rd}
+%\usepackage{graphicx}
+\usepackage[dvipdfm]{graphicx} %for pdf output
+\usepackage{color}
+\usepackage{amsmath,amssymb}
+\usepackage{fancybox} % this is needed for \VerbatimEnvironment
+\usepackage{listings}
+\definecolor{mygray}{rgb}{0.96,0.96,0.96}
+\definecolor{myblack}{rgb}{0,0,0}
+\lstset{language=R, captionpos=b}
+\lstset{rulesepcolor=\color{myblack}, backgroundcolor=\color{mygray}}
+\lstset{showspaces=false, showtabs=false, showstringspaces=false}
+\lstset{basicstyle=\ttfamily\scriptsize}
+
+
+\setlength{\topmargin}{-3cm}
+\textwidth=17.5cm
+\textheight=26cm
+\oddsidemargin=-1cm
+
+\def\de{{\rm d}}
+\def\diag{\mathop{\rm diag}\nolimits}
+\begin{document}
+
+
+\maketitle
+
+%Abst
+
+\SweaveInput{nakamacro-aism20}
+\SweaveInput{yuima_1}
+\SweaveInput{yuima_2}
+\SweaveInput{yuima_3}
+\SweaveInput{yuima_4}
+\SweaveInput{yuima_5}
+\SweaveInput{bib}
+
+
+\end{document}
Added: pkg/yuimadocs/inst/doc/yuima_1.Rnw
===================================================================
--- pkg/yuimadocs/inst/doc/yuima_1.Rnw (rev 0)
+++ pkg/yuimadocs/inst/doc/yuima_1.Rnw 2009-11-07 18:50:57 UTC (rev 11)
@@ -0,0 +1,53 @@
+\section{Introduction}
+
+\def\yuima{{\bf yuima}}
+
+The package $\yuima$ provides
+an object-oriented programming environment
+for simulation and statistical inference
+for stochastic processes by R.
+On this framework,
+$\yuima$ also supplies various functions
+to execute simulation and statistical analysis.
+Both categories of procedures may depend 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.
+%It is therefore natural to
+$\yuima$ 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.
+For example, the asymptotic expansion method for computing
+option prices (i.e., expectation of an irregular functional of
+a stochastic process) provides precise approximation values
+instantaneously, taking advantage of the analytic approach,
+but the formula has a long expression like more than one page!
+$\yuima$ delivers up-to-date methods as a package
+onto the desk of the user working
+with simulation and/or statistics for stochastic differential equations.
+
+
+Sampled data from a continuous-time process features
+the time stamps as well as the positions of the object.
+It is requiring a new theory of estimation.
+$\yuima$ framework can apply multi-dimensional time stamps
+of tick data and provides diverse functions handling such kind data
+to support statistical analysis.
+
Added: pkg/yuimadocs/inst/doc/yuima_2.Rnw
===================================================================
--- pkg/yuimadocs/inst/doc/yuima_2.Rnw (rev 0)
+++ pkg/yuimadocs/inst/doc/yuima_2.Rnw 2009-11-07 18:50:57 UTC (rev 11)
@@ -0,0 +1,218 @@
+
+\section{Yuima object}
+The main object in our framework is the \code{yuima} object.
+The \code{yuima} object contains slots named \code{model}, \code{data},
+\code{sampling}, and \code{characteristic}.
+
+\begin{center}
+\begin{tabular}{|c|}
+\hline
+\colorbox{black}{\textcolor{white}{ \code{yuima} \bf object}}\\
+\hline
+\hline
+\code{model}\\
+\hline
+\code{data}\\
+\hline
+\code{sampling}\\
+\hline
+\code{characteristic}\\
+\hline
+\end{tabular}
+\end{center}
+
+
+The slot \code{data} will contain
+a time series data and the slot \code{sampling} contains a description of
+sampling methods and/or parameters of the data collection(e.g. poisson sampling with rate $p=0.1$.).
+
+The \code{yuima.model} object is intended to be a description of model
+independent of the actual data.
+For example, in model selection or hypothesis testing for SDEs,
+it is possible to have different models to fit to the same set of data.
+
+On the contrary, for the same model specified in \code{yuima.model},
+different time series can be fitted.
+In this case the vector \code{model} has length 1, and
+the vector of \code{data} has length greater than 1.
+
+The \code{yuima.characteristic} object is prepared for specifying the characteristics of the \code{yuima} object such as
+the number of equations in the model and data, time and space unit scale for simulation or data.
+For now, this slot is mainly for keeping consistency of the other three slots.
+%For example, the number of pathes specified in the \code{yuima.model} object and the actual data in \code{yuima.data} is
+%not the same,
+
+\code{setYuima} is a constructor for \code{yuima} class,
+which does not initialize but only creates four NULL slots:
+
+\begin{itemize}
+\item{data:}{ an object of class \code{yuima.data}}
+\item{model:}{ an object of class \code{yuima.model}}
+\item{sampling:}{ an object of class \code{yuima.sampling}}
+\item{characteristic:}{ an object of class \code{yuima.characteristic}}
+\end{itemize}
+%We show how to construct \code{yuima} object in the next section.
+
+%%Note that \code{yuima.sampling} is not implemented yet.
+
+%-----------------------------------
+ \begin{figure}[!htb]
+ \centering
+ \includegraphics[width=10cm,clip]{./EPS/functionDefinition.eps}
+ \caption{function relation}
+ \end{figure}
+ %-----------------------------------
+%-----------------------------------
+ \begin{figure}[!htb]
+ \centering
+ \includegraphics[width=5cm,clip]{./EPS/classDefinition.eps}
+ \caption{class definition}
+ \end{figure}
+ %-----------------------------------
+
+\subsection{yuima.model}
+\code{yuima.model} specifies the SDE of interest. The constructor \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 \code{yuima.model} contains several slots listed below.
+To see inside its structure, we use the \R{} command \code{str}.
+\begin{itemize}
+\item \code{drift} is an \R{} expression which contains the drift specification.
+\item \code{diffusion} is itself a list of 1 slot which describes the diffusion
+ coefficient relative to first noise.
+
+\item \code{parameter} which is a short name for ``parameters'' which
+ is a list of objects.
+
+\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{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{state.variable} and \code{time.variable}, by default,
+ are assumed to be $x$ and $t$ but the user can freely choose them.
+ The \code{yuima.model} function assumes that the user either use 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{noise.number} indicates the number of sources of noise.
+
+\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.
+
+\end{itemize}
+We show an example to describe the following simple stochastic differential equation
+
+\begin{equation}
+\mbox{d}X_t=3X_t\mbox{d}t+\mbox{d}W_t.
+\end{equation}
+
+<<print=FALSE,echo=TRUE>>=
+library(yuima)
+@
+
+<<echo=TRUE>>=
+mod <- setModel(drift="3*x")
+str(mod)
+@
+
+\subsection{yuima.data}
+\code{yuima.data} is a class to represent the data to be used. The constructor \code{setData} is
+used to give an observed data to the \code{yuima.data} object. \code{get.zoo.data} returns the \code{zoo.data} field of
+an \code{yuima.data} object.
+As usual class methods, \code{yuima.data} provides \code{length}, \code{plot} and
+\code{dim} functions to get necessary information.
+
+ The \code{zoo.data} slot may contain different kind of observed data.
+ In the \R{} system, there are different kinds of objects to describe time series data.
+ Here is a very short list of possibilities
+ \begin{center}
+ \begin{tabular}{r l}
+ \code{ts}, \code{mts}: & regularly spaced time series 1-dim or multidim\\
+ \code{its}: & irregularly spaced time series\\
+ \code{tseries}: & time index is real time of financial markets \\
+ \code{zoo}: & more flexibility in the management of the index of the time series.\\
+ \end{tabular}
+ \end{center}
+ Because \R{} is used in many contexts, including finance, and users adopt one of
+ the above times series objects because they are used to it,
+ our framework accepts any of the above as observed data and treat the inference accordingly.
+ \\
+ \\
+ An object of \code{yuima.data} contains several slots listed below.
+ To see inside its structure we use the \R{} command \code{str}.
+ \begin{itemize}
+ \item \code{original.data} is a copy of the original data. (it can be any type of data.)
+ \item \code{zoo.data} is an bundled \code{zoo} version object of class \code{yuima.data} or \code{yuima}.
+ \end{itemize}
+
+ We show an example of \code{yuima.data} object.
+
+<<print=FALSE,echo=TRUE>>=
+library(zoo)
+@
+
+<<fig=TRUE,echo=TRUE,width=7,height=7,results=hide>>=
+X <- ts(matrix(rnorm(200),100,2))
+mydata <- setData(X)
+str(get.zoo.data(mydata))
+#plot(get.zoo.data(mydata))
+dim(mydata)
+length(mydata)
+plot(mydata)
+@
+
+Then, we show how to construct \code{yuima} object with \code{yuima.data} and \code{yuima.model} object.
+
+We first construct a two-dimensional homogeneous diffusion process:
+\begin{eqnarray*}
+\mbox{d} X_{1t} &=&
+- 3 X_{1t} \mbox{d}t + \mbox{d} W_{1t}, \\
+\mbox{d}X_{2t} &=&
+-(X_{1t} + 2 X_{2t}) \mbox{d}t + 0.5 \mbox{d} W_{1t} + \mbox{d}W_{2t}.
+\end{eqnarray*}
+<<echo=TRUE>>=
+# we set the drift coefficient
+a <- c("-3*x1","-x1-2*x2")
+# and also the diffusion coefficient
+b <- matrix(c("1","0.5","0","1"),2,2)
+# Then set
+mod2 <- setModel(drift = a, diffusion = b, solve.variable = c("x1","x2"))
+str(mod2)
+@
+
+Then, we can combine the \code{yuima.model} and \code{yuima.data}
+objects to a \code{yuima} object.
+<<echo=TRUE>>=
+yuima.obj <- setYuima(data=mydata,model=mod2)
+str(yuima.obj)
+@
+
+\subsection{Realworld data}
+In this subsection, we show a simple example of handling the realworld data in
+the yuima framework.
+<<fig=TRUE,echo=TRUE,width=7,height=7>>=
+## Use TTR package to get real market data from Yahoo! Finance
+require(TTR)
+sbux<-getYahooData("SBUX", start=20080401, end=20090331)
+cbou<-getYahooData("CBOU", start=20080401, end=20090331)
+sd <- setData(list(sbux,cbou))
+real.yuima <- setYuima(data=sd)
+summary(real.yuima)
+head(get.zoo.data(real.yuima)[[1]])
+head(get.zoo.data(real.yuima)[[2]])
+
+plot(get.zoo.data(real.yuima)[[1]]$High,main="SBUX_High")
+@
Added: pkg/yuimadocs/inst/doc/yuima_3.Rnw
===================================================================
--- pkg/yuimadocs/inst/doc/yuima_3.Rnw (rev 0)
+++ pkg/yuimadocs/inst/doc/yuima_3.Rnw 2009-11-07 18:50:57 UTC (rev 11)
@@ -0,0 +1,138 @@
+\section{Simulation}
+\subsection{Diffusion Process}
+\code{simulate} is a function to solve SDE using the Euler-Maruyama method.
+Here we illustrate how to generate a simulated path of the solution of
+a given SDE. Let us consider a 3-dim SDE
+
+ \begin{equation} \label{sabr}
+ \begin{cases}
+ & \mbox{d} X_t^1 = X_t^2 \left| X_t^1 \right|^{2/3} \mbox{d}W_t^1, \\
+ & \mbox{d}X_t^2 = g(t)\mbox{d}X_t^3, \\
+ & \mbox{d}X_t^3 = X_t^3( \mu \mbox{d}t + \sigma (\rho \mbox{d}W_t^1 + \sqrt{1-\rho^2}
+ \mbox{d}W_t^2))
+ \end{cases},
+ \ \ \ \ (X_0^1, X_0^2, X_0^3) = (1.0,0.1,1.0)
+ \end{equation}
+ with $\mu = 0.1, \sigma = 0.2, \rho = -0.7 $ and $g(t) = 0.4 + (0.1 + 0.2t) e^{-2t}$, for
+ example, where $W = (W^1, W^2)$ is a 2-dim standard Brownian motion.
+ The solution $X = (X^1, X^2, X^3)$ takes values on $\mathbb{R}_+^3$.
+ This is a stochastic integral equation defined as
+ \begin{equation}
+ X_t = X_0 + \int_0^t b(s,X_s) \mbox{d}s + \int_0^t c(s,X_s) \mbox{d}W_s
+ \end{equation}
+ with
+ \[
+ b(s,x) = \left( \begin{array}{@{\,}c@{\,}} 0 \\ g(s) \mu x^3 \\ \mu x^3 \end{array} \right),
+ \ \
+ c(s,x) = \left( \begin{array}{@{\,}cc@{\,}} x^2 | x^1|^{2/3} & 0 \\
+ g(s)\sigma \rho x^3 & g(s)\sigma \sqrt{1-\rho^2} x^3 \\
+ \sigma \rho x^3 & \sigma \sqrt{1-\rho^2} x^3
+ \end{array} \right)
+ \]
+ for $ x = (x^1, x^2, x^3)$.
+ The Euler-Maruyama method with discretization scheme
+ \begin{equation}
+ 0 = \tau_0 < \tau_1 < \cdots < \tau_j < \tau_{j+1} < \cdots
+ \end{equation}
+ refers to an approximative solution given by
+ \begin{equation}
+ \tilde{X}_{\tau_{j+1}} = \tilde{X}_{\tau_j} +
+ b(\tau_j,\tilde{X}_{\tau_j})(\tau_{j+1} - \tau_j) +
+ c(\tau_j, \tilde{X}_{\tau_j} )(W_{\tau_{j+1}}- W_{\tau_j}), \ \ j \geq 0.
+ \end{equation}
+ If interval $\tau_{j+1} - \tau_j$ is independent to $\{W_t\}_{t \geq \tau_j}$ for each $j$,
+ $W_{\tau_{j+1}} - W_{\tau_j}$ has the same distribution as $\sqrt{\tau_{j+1} - \tau_j}N_j$,
+ where $N_j$ is an iid sequence of standard normal variables.
+ The usual discretization scheme is an equidistant one such as
+ \begin{equation}
+ \tau_{j+1} - \tau_j = \text{Terminal}/\text{division}
+ \end{equation}
+ for all $j \geq 0$.
+
+<<fig=TRUE,echo=TRUE,width=7,height=7>>=
+ set.seed(123)
+
+ mu <- 0.1
+ sig <- 0.2
+ rho <- -0.7
+
+ g <- function(t){
+ return(0.4 + (0.1 + 0.2*t)* exp(-2*t))
+ }
+
+ f1 <- function(t, x1, x2, x3) {
+ ret <- 0
+ if(x1 > 0 && x2 > 0) ret <- x2*exp(log(x1)*2/3)
+ return(ret)
+ }
+
+ f2 <- function(t, x1, x2, x3) {
+ ret <- 0
+ if(x3 > 0) ret <- rho*sig*x3
+ return(ret)
+ }
+
+ f3 <- function(t, x1, x2, x3) {
+ ret <- 0
+ if(x3 > 0) ret <- sqrt(1-rho^2)*sig*x3
+ return(ret)
+ }
+
+ diff.coef.matrix <- matrix(c("f1(t,x1,x2,x3)", "f2(t,x1,x2,x3) * g(t)",
+ "f2(t,x1,x2,x3)", "0", "f3(t,x1,x2,x3)*g(t)", "f3(t,x1,x2,x3)"), 3, 2)
+
+ sabr.mod <- setModel(drift = c("0", "mu*g(t)*x3", "mu*x3"), diffusion = diff.coef.matrix,
+ state.variable = c("x1", "x2", "x3"), solve.variable = c("x1", "x2", "x3"))
+
+ # Change the parameters below if needed
+ yuima.samp <- setSampling(Terminal=1,division=1000)
+
+ # make yuima object
+ yuima <- setYuima(model=sabr.mod, sampling=yuima.samp)
+
+ # Solve SDEs using Euler-Maruyama method
+ yuima <- simulate(yuima, xinit=c(1.0,0.1,1.0))
+ plot(yuima)
+@
+
+ \subsection{Space-discretizing scheme}
+
+ \code{simulate} provides also space-discretized Euler-Maruyama method to solve SDE.
+ This is at the moment designed for 1 factor SDE, i.e., the case the driving
+ Brownian motion $W$ is 1-dimensional. Here discretization scheme $\{\tau_j\}$
+ is defined as
+
+\begin{equation}
+\tau_0 = 0, \ \ \tau_{j+1} = \inf\{ t > \tau_j; |W_t - W_{\tau_j}| = \epsilon \}
+\end{equation}
+for each $j \geq 0$. \code{simulate} takes
+\begin{equation}
+\epsilon^2 = \text{Terminal} / \text{division},
+\end{equation}
+which coinsides with the mean of the interval $\tau_{j+1} - \tau_j$.
+This space-discretizing scheme is 3 times efficient than the usual time-discretizing
+one in the sense of mean-squared error.
+
+Let us consider the previous model (\ref{sabr}) with $\rho = -1$.
+
+<<fig=TRUE,echo=TRUE,width=7,height=7>>=
+
+set.seed(123)
+
+rho <- -1
+
+diff.coef.matrix <- c("f1(t,x1,x2,x3)", "f2(t,x1,x2,x3) * g(t)",
+ "f2(t,x1,x2,x3)")
+
+sabr.mod <- setModel(drift = c("0", "mu*g(t)*x3", "mu*x3"), diffusion = diff.coef.matrix,
+ state.variable = c("x1", "x2", "x3"), solve.variable = c("x1", "x2", "x3"))
+
+yuima <- setYuima(model=sabr.mod, sampling=yuima.samp)
+
+#yuima <- simulate(yuima)
+yuima.sd <- simulate(yuima, space.discretized = TRUE, xinit=c(1.0,0.1,1.0))
+
+#plot(yuima)
+plot(yuima.sd)
+@
+
Added: pkg/yuimadocs/inst/doc/yuima_4.Rnw
===================================================================
--- pkg/yuimadocs/inst/doc/yuima_4.Rnw (rev 0)
+++ pkg/yuimadocs/inst/doc/yuima_4.Rnw 2009-11-07 18:50:57 UTC (rev 11)
@@ -0,0 +1,300 @@
+\section{cce}
+ cce is a function to estimate the covariance between two It\^o
+ processes when they are observed at discrete times nonsynchronously.
+ It can also apply to irregularly sampled one-dimensional data as a special case.
+
+\subsection{Nonsynchronous covariance estimator}
+Suppose that two It\^o processes are observed
+only at discrete times in a nonsynchronous manner.
+We are interested in estimating the covariance of the two processes accurately
+in such a situation.
+This type of problem arises typically in high-frequency financial time series.
+
+Let $T \in (0,\infty)$ be a terminal time for possible observations.
+We consider a two dimensional It\^o process $(X^1,X^2)$
+satisfying the stochastic differential equations
+\begin{eqnarray*}
+\mbox{d}X_t^l &=&\mu_t^l \mbox{d}t + \sigma_t^l \mbox{d}W_t^l ,\quad
+t\in[0,T]\\
+X_0^l &=& x_0^l \notag
+\end{eqnarray*}
+for $l=1,2$.
+Here $W^l$ denote standard Wiener processes with a
+progressively measurable correlation process
+$\mbox{d}\langle W_1 , W_2 \rangle_t = \rho_t \mbox{dt}$,
+$\mu^l_t$ and $\sigma_t^l$ are progressively measurable processes,
+and $x^l_0$ are initial random variables independent of $(W^1,W^2)$.
+Diffusion type processes are in the scope but this model
+can express more sophisticated stochastic structures.
+
+
+The process $X^l$ is supposed to be observed at
+over the increasing sequence of times
+$T^{l,i}$ $(i\in\mathbb{Z}_{\geq0})$ starting at $0$, up to time T.
+Thus, the observables are $(T^{l,i},X^{l,i})$ with $T^{l,i}\leq T$.
+Each $T^{l,j}$ may be a stopping time, so possibly depends on
+the history of $(X^1,X^2)$ as well as the precedent stopping times.
+Two sequences of stopping times
+$T^{1,i}$ and $T^{2,j}$ are {\it nonsynchronous}, and
+irregularly spaced, in general.
+In particular, cce can apply to estimation of the quadratic variation of a single
+stochastic process sampled regularly/irregularly.
+
+
+The parameter of interest is the quadratic covariation between $X^1$ and $X^2$:
+
+\begin{equation}
+\theta=
+\langle X^1 , X^2 \rangle_T = \int_0^T \sigma_t^1 \sigma_t^2 \rho_t \mbox{d}t.
+\end{equation}
+The target variable $\theta$ is random in general.
+
+It can be estimated with the nonsynchronous covariance estimator
+(Hayashi-Yoshida estimator)
+\begin{equation}
+U_n= \sum_{i,j:T^{1,i}\leq T, T^{2,j}\leq T} (X_{T^{1,i}}^1-X_{T^{1,i-1}}^1)(X_{T^{2,j}}^2-X_{T^{2,j-1}}^2)
+1_{\{ (T^{1,i-1},T^{1,i}]
+\bigcap (T^{2,j-1},T^{2,j}] \neq \emptyset \}}.
+\end{equation}
+That is, the product of any pair of increments $(X_{T^{1,i}}^1-X_{T^{1,i-1}}^1)$ and
+$(X_{T^{2,j}}^2-X_{T^{2,j-1}}^2)$ will make a contribution
+to the sum only when
+the respective observation intervals $(T^{1,i-1},T^{1,i}] $ and $ (T^{2,j-1},T^{2,j}] $ are overlapping
+with each other.
+%The estimator $U_n$ was proposed by Hayashi and Yoshida and
+It is known that $U_n$ is consist and has asymptotically mixed normal
+distribution
+as $n\to\infty$ if the maximum length between two consecutive observing times
+tends to $0$.
+See \cite{hay-yos05,hay-yos04,hay-yos06,hay-yos08} for details.
+
+\subsection{Example: data generation and estimation by yuima package}
+We will demonstrate how to apply cce function to
+nonsynchronous high-frequency data by simulation.
+As an example, consider a two dimensional stochastic process
+$(X^1_t,X^2_t)$ satisfying the stochastic differential equation
+\begin{equation}
+\begin{split}
+\mbox{d}X^1_t = \sigma_{1,t} \mbox{d}B^1_t, \\
+\mbox{d}X^2_t = \sigma_{2,t} \mbox{d}B^2_t.
+\end{split}
+\end{equation}
+Here $B^1_t$ and $B^2_t$ denote two standard Wiener processes,
+however they are correlated as
+\begin{eqnarray}
+B^1_t &=& W^1_t,\\
+B^2_t &=& \int_0^t \rho_s \mbox{d} W^1_s +
+ \int_0^t \sqrt{ 1- \rho_s^2} \mbox{d} W^2_s,
+\end{eqnarray}
+where $W^1_t$ and $W^2_t$ are independent Wiener processes, and
+$\rho_t$ is the correlation function between $B^1_t$ and $B^2_t$.
+We consider $\sigma_{l,t},l=1,2$ and $\rho_t$
+of the following form in this example:
+
+\begin{eqnarray*}
+ \sigma_{1,t} &=& \sqrt{1+t}, \\
+ \sigma_{2,t} &=& \sqrt{1+t^2}, \\
+ \rho_t &=& \frac{1}{\sqrt{2}}.
+\end{eqnarray*}
+
+To simulate the stochastic process $(X^1_t,X^2_t)$, we first build the model
+by setModel as follows.
+It should be noted that the method of generating
+nonsynchronous data can be replaced by a simpler one but
+we will take a general approach here to demonstrate
+a usage of
+the yuima comprehensive package for simulation and estimation of
+stochastic processes.
+
+<<print=FALSE,echo=TRUE>>=
+# diffusion coefficient for process 1
+diff.coef.1 <- function(t,x1=0, x2=0) sqrt(1+t)
+# diffusion coefficient for process 2
+diff.coef.2 <- function(t,x1=0, x2=0) sqrt(1+t^2)
+# correlation
+cor.rho <- function(t,x1=0, x2=0) sqrt(1/2)
+# coefficient matrix for diffusion term
+diff.coef.matrix <- matrix( c( "diff.coef.1(t,x1,x2)", "diff.coef.2(t,x1,x2) * cor.rho(t,x1,x2)", "", "diff.coef.2(t,x1,x2) * sqrt(1-cor.rho(t,x1,x2)^2)"),2,2)
+# Model sde using yuima.model
+cor.mod <- setModel(drift = c("",""), diffusion = diff.coef.matrix, solve.variable=c("x1","x2"))
+@
+
+The parameter we want to estimate is the quadratic covariation between $X_1$ and $X_2$:
+
+\begin{equation}
+\theta = \langle X_1, X_2 \rangle_T =
+ \int_0^T \sigma_{1,t} \sigma_{2,t} \rho_t \mbox{d} t.
+\end{equation}
+Later, we will compare estimated values with
+the true value of $\theta$ given by
+<<echo=TRUE>>=
+CC.theta <- function( T, sigma1, sigma2, rho)
+{
+ tmp <- function(t) return( sigma1(t) * sigma2(t) * rho(t) )
+ integrate(tmp,0,T)
+}
+# Cumulative Covariance
+theta <- CC.theta(T=1, sigma1=diff.coef.1, sigma2=diff.coef.2, rho=cor.rho)$value
+cat(sprintf("theta=%5.3f\n",theta))
+@
+For the sampling scheme, we will consider the independent oisson sampling.
+That is, each configuration of the sampling times $T^{l,i}$ is realized
+as the Poisson random measure with intensity $np_l$,
+and the two random measures are independent each other as well as
+the stochastic processes.
+Then it is known that
+
+
+\begin{equation}
+n^{1/2} ( U_n -\theta) \rightarrow N(0,c),
+\end{equation}
+as $n\to\infty$,
+where
+\begin{equation}
+ c = \left( \frac{2}{p_1} + \frac{2}{p_2} \right)
+\int_0^T \left( \sigma_{1,t} \sigma_{2,t} \right)^2 \mbox{d}t +
+\left(
+\frac{2}{p_1} + \frac{2}{p_2} - \frac{2}{p_1 + p_2}
+\right)
+\int_0^T
+\left(
+\sigma_{1,t} \sigma_{2,t} \rho_t
+\right)^2 \mbox{d} t.
+\end{equation}
+
+We use a function \code{poisson.random.sampling} to get observations
+by the Poisson sampling.
+This function returns a matrix combined times and observations of
+Poisson random sampling.
+
+<<echo=TRUE>>=
+set.seed(111)
+
+Terminal <- 1
+p1 <- 0.2
+p2 <- 0.3
+rate<-c(p1, p2)
+n <- 1000
+# Change to n<-1000
+division <- 1200
+# Change too division <- 1200
+
+yuima.samp <- setSampling(Terminal=Terminal,division=division)
+yuima <- setYuima(model=cor.mod, sampling=yuima.samp)
+
+# solve SDEs using Euler-Maruyama method
+Data <- simulate(yuima)
+
+sampledData <- poisson.random.sampling(Data,rate=rate, n=n)
+@
+\code{cce} takes the sample and returns
+an estimate of the quadratic covariation.
+<<cce, echo=TRUE>>=
+cce(sampledData)[1,2] # CumulativeCovarianceEstimator
+@
+and in case of non missing data
+<<cce2, echo=TRUE>>=
+cce(Data)[1,2] # CumulativeCovarianceEstimator
+@
+The last value is just for reference; it may bear another kind error
+due to approximation error to the Poisson sampling.
+
+
+In what follows, we will evaluate the difference between
+the estimated value and the true value, as well as the asymptotic variance,
+and compare the histogram of the estimation error after scaling
+with the limiting Gaussian distribution
+(Figure \ref{fig:histogram}).
+The function \texttt{var.c} calculates the theoretical
+asymptotic variance of the scaled estimation error.
+
+<<cctheta2, echo=TRUE>>=
+#error <- NULL # difference between true parameter and estimator
+
+
+# Cumulative Covariance
+theta <- CC.theta(T=Terminal, diff.coef.1, diff.coef.2, cor.rho)$value
+cat(sprintf("The true value of theta=%5.3f\n",theta))
+
+
+
+# asymptotic variance
+var.c <- function(T, p1,p2, sigma1, sigma2, rho)
+{
+ tmp_integrand1 <- function(t) (sigma1(t) * sigma2(t))^2
+ i1 <- integrate(tmp_integrand1,0,T)
+ tmp_integrand2 <- function(t) (sigma1(t) * sigma2(t) * rho(t))^2
+ i2 <- integrate(tmp_integrand2,0,T)
+ 2*(1/p1 + 1/p2)* i1$value + 2*(1/p1+1/p2 - 1/(p1+p2)) * i2$value
+}
+
+vc <- var.c(T=Terminal, p1, p2, diff.coef.1, diff.coef.2, cor.rho)
+sd <- sqrt(vc/n)
+cat(sprintf("The standard deviation by normal approximation =%5.3f\n",sd))
+
+
+
+# simulation
+N <- 100
+U <- NULL
+for(i in 1:N)
+{
+ # solve SDEs using Euler-Maruyama method
+ yuima.tmp <- simulate(yuima)
+ yuimaData <- poisson.random.sampling(yuima.tmp,rate=rate,n=n)
+ Un <- cce(yuimaData)[1,2] # CumulativeCovarianceEstimator
+ U <- append(U, Un)
+}
+
+cat(sprintf("The sample mean =%5.3f\n", mean(U)))
+cat(sprintf("The sample standard deviation =%5.3f\n", sqrt(var(U))))
+@
+
+
+\begin{figure}[h]
+\centering
+<<fig1, fig=TRUE,echo=TRUE,width=10,height=10,results=hide>>=
+# plot sampled data histogram and Gaussian distribution (theoretical result).
+library(KernSmooth)
+h <- dpih(U)
+bins <- seq(min(U)-5, max(U)+5+h,by=h)
+xrange <- c(min(U)-5-h, max(U)+5+h)
+yrange<-c(0,max(hist(U)$density))
+hist(U,breaks=bins,xlab="",ylab="",xlim=xrange,ylim=yrange,probability=T)
+curve(dnorm(x,mean=theta,sd=sd),ylab="",xlim=xrange,ylim=yrange,col="red", add=TRUE)
+@
+\caption{A sampled data histogram and Gaussian distribution(theoretical result)}
+\label{fig:histogram}
+\end{figure}
+
+
+
+
+\begin{thebibliography}{99}
+\bibitem{hay-yos05}
+Hayashi, T., Yoshida, N.:
+On Covariance Estimation of Non-synchronously Observed Diffusion Processes.
+Bernoulli, {\bf 11}, 2, 359-379 (2005)
+
+
+\bibitem{hay-yos04}
+Hayashi, T., Yoshida, N.:
+Asymptotic normality of a covariance estimator for nonsynchronously observed diffusion processes.
+Annals of the Institute of Statistical Mathematics,
+{\bf 60}, 2, %/ 2008”N6ŒŽ
+367-406 (2008)
+
+\bibitem{hay-yos06}
+Hayashi, T., Yoshida, N.:
+Nonsynchronous Covariance Estimator and Limit Theorem.
+Research Memorandum No.1020, Institute of Statistical Mathematics, 2006
+
+\bibitem{hay-yos08}
+Hayashi, T., Yoshida, N.:
+Nonsynchronous Covariance Estimator and Limit Theorem II.
+Research Memorandum No.1067, Institute of Statistical Mathematics, 2008
+
+\end{thebibliography}
+
+
+
More information about the Yuima-commits
mailing list