[Yuima-commits] r206 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Oct 5 06:44:45 CEST 2012
Author: kyuta
Date: 2012-10-05 06:44:45 +0200 (Fri, 05 Oct 2012)
New Revision: 206
Added:
pkg/yuima/R/bns.test.R
pkg/yuima/R/mpv.R
pkg/yuima/man/bns.test.Rd
pkg/yuima/man/mpv.Rd
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NEWS
Log:
add mpv.R, bns.test.R, mpv.Rd, bns.test.Rd
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2012-07-26 07:59:00 UTC (rev 205)
+++ pkg/yuima/DESCRIPTION 2012-10-05 04:44:45 UTC (rev 206)
@@ -1,8 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project package (unstable version)
-Version: 0.1.196
-Date: 2012-05-24
+Version: 0.1.197Date: 2012-10-05
Depends: methods, zoo, stats4, utils
Suggests: cubature, mvtnorm
Author: YUIMA Project Team.
Modified: pkg/yuima/NEWS
===================================================================
--- pkg/yuima/NEWS 2012-07-26 07:59:00 UTC (rev 205)
+++ pkg/yuima/NEWS 2012-10-05 04:44:45 UTC (rev 206)
@@ -0,0 +1 @@
+2012/10/05: add mpv.R, bns.test.R, mpv.Rd, bns.test.Rd
\ No newline at end of file
Added: pkg/yuima/R/bns.test.R
===================================================================
--- pkg/yuima/R/bns.test.R (rev 0)
+++ pkg/yuima/R/bns.test.R 2012-10-05 04:44:45 UTC (rev 206)
@@ -0,0 +1,83 @@
+
+##########################################################################
+# Barndorff-Nielsen and Shephard jump test
+##########################################################################
+
+setGeneric("bns.test",
+ function(yuima,r=rep(1,4),type="standard",adj=TRUE)
+ standardGeneric("bns.test"))
+
+setMethod("bns.test",signature(yuima="yuima"),
+ function(yuima,r=rep(1,4),type="standard",adj=TRUE)
+ bns.test(yuima at data,r,type,adj))
+
+setMethod("bns.test",signature(yuima="yuima.data"),
+ function(yuima,r=rep(1,4),type="standard",adj=TRUE){
+
+ # functions for computing the test statistics
+
+ bns.stat_standard <- function(yuima,r=rep(1,4)){
+
+ JV <- mpv(yuima,2)-mpv(yuima)
+ IQ <- mpv(yuima,r)
+
+ return(JV/sqrt((pi^2/4+pi-5)*IQ))
+ }
+
+ bns.stat_ratio <- function(yuima,r=rep(1,4),adj=TRUE){
+
+ bpv <- mpv(yuima)
+ RJ <- 1-bpv/mpv(yuima,2)
+ avar <- mpv(yuima,r)/bpv^2
+
+ if(adj){
+ avar[avar<1] <- 1
+ }
+
+ return(RJ/sqrt((pi^2/4+pi-5)*avar))
+ }
+
+ bns.stat_log <- function(yuima,r=rep(1,4),adj=TRUE){
+
+ bpv <- mpv(yuima)
+ RJ <- log(mpv(yuima,2)/bpv)
+ avar <- mpv(yuima,r)/bpv^2
+
+ if(adj){
+ avar[avar<1] <- 1
+ }
+
+ return(RJ/sqrt((pi^2/4+pi-5)*avar))
+ }
+
+ # main body of the test procedure
+
+ data <- get.zoo.data(yuima)
+ d.size <- length(data)
+
+ n <- integer(d.size)
+ for(d in 1:d.size){
+ n[d] <- sum(!is.na(as.numeric(data[[d]])))
+ if(n[d]<2) {
+ stop("length of data (w/o NA) must be more than 1")
+ }
+ }
+
+ switch(type,
+ "standard"="<-"(bns,bns.stat_standard(yuima,r)),
+ "ratio"="<-"(bns,bns.stat_ratio(yuima,r,adj)),
+ "log"="<-"(bns,bns.stat_log(yuima,r,adj)))
+
+ bns <- sqrt(n)*bns
+
+ result <- vector(d.size,mode="list")
+ for(d in 1:d.size){
+ p <- pnorm(bns[d],lower.tail=FALSE)
+ result[[d]] <- list(statistic=c(BNS=bns[d]),p.value=p,
+ method="Barndorff-Nielsen and Shephard jump test",
+ data.names=paste("x",d,sep=""))
+ class(result[[d]]) <- "htest"
+ }
+
+ return(result)
+ })
Added: pkg/yuima/R/mpv.R
===================================================================
--- pkg/yuima/R/mpv.R (rev 0)
+++ pkg/yuima/R/mpv.R 2012-10-05 04:44:45 UTC (rev 206)
@@ -0,0 +1,55 @@
+
+#########################################################################
+# Realized multipower variation
+#########################################################################
+
+setGeneric("mpv",function(yuima,r=c(1,1))standardGeneric("mpv"))
+
+setMethod("mpv",signature(yuima="yuima"),function(yuima,r=c(1,1))mpv(yuima at data,r))
+
+setMethod("mpv",signature(yuima="yuima.data"),function(yuima,r=c(1,1)){
+
+ data <- get.zoo.data(yuima)
+ d.size <- length(data)
+
+ result <- double(d.size)
+
+ if(is.numeric(r)){
+ for(d in 1:d.size){
+
+ X <- as.numeric(data[[d]])
+ idt <- which(is.na(X))
+ if(length(idt>0)){
+ X <- X[-idt]
+ }
+ if(length(X)<2) {
+ stop("length of data (w/o NA) must be more than 1")
+ }
+
+ abs.diffX <- abs(diff(X))
+ tmp <- rollapplyr(abs.diffX,length(r),FUN=function(x)prod(x^r))
+ result[d] <- length(abs.diffX)^(sum(r)/2-1)*sum(tmp)/
+ prod(2^(r/2)*gamma((r+1)/2)/gamma(1/2))
+
+ }
+ }else{
+ for(d in 1:d.size){
+
+ X <- as.numeric(data[[d]])
+ idt <- which(is.na(X))
+ if(length(idt>0)){
+ X <- X[-idt]
+ }
+ if(length(X)<2) {
+ stop("length of data (w/o NA) must be more than 1")
+ }
+
+ abs.diffX <- abs(diff(X))
+ tmp <- rollapplyr(abs.diffX,length(r[[d]]),FUN=function(x)prod(x^r[[d]]))
+ result[d] <- length(abs.diffX)^(sum(r[[d]])/2-1)*sum(tmp)/
+ prod(2^(r[[d]]/2)*gamma((r[[d]]+1)/2)/gamma(1/2))
+ }
+ }
+
+ return(result)
+})
\ No newline at end of file
Added: pkg/yuima/man/bns.test.Rd
===================================================================
--- pkg/yuima/man/bns.test.Rd (rev 0)
+++ pkg/yuima/man/bns.test.Rd 2012-10-05 04:44:45 UTC (rev 206)
@@ -0,0 +1,118 @@
+\name{bns.test}
+\alias{bns.test,list-method}
+%- Also NEED an '\alias' for EACH other topic documented here.
+\title{
+%% ~~function to do ... ~~
+Barndorff-Nielsen and Shephard's Test for the Presence of Jumps Using Bipower Variation
+}
+\description{
+%% ~~ A concise (1-5 lines) description of what the function does. ~~
+Tests the presence of jumps using the statistic proposed in Barndorff-Nielsen and Shephard (2004,2006) for each components.
+}
+\usage{
+bns.test(yuima, r = rep(1, 4), type = c("standard","log","ratio"), adj = TRUE)
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+ \item{yuima}{
+%% ~~Describe \code{yuima} here~~
+an object of \code{\link{yuima-class}} or \code{\link{yuima.data-class}}.
+}
+ \item{r}{
+%% ~~Describe \code{r} here~~
+a vector of non-negative numbers or a list of vectors of non-negative numbers.
+}
+ \item{type}{
+%% ~~Describe \code{type} here~~
+type of test statistics to use. \code{standard} is default.
+}
+ \item{adj}{
+%% ~~Describe \code{adj} here~~
+logical; if \code{TRUE}, the muximum adjustment suggested in Barndorff-Nielsen and Shephard (2004) is applied to the test statistic when \code{type} is equal to either \dQuote{\code{log}} or \dQuote{\code{ratio}}.
+}
+}
+\details{
+%% ~~ If necessary, more details than the description above ~~
+For the \code{i}-th component, the test statistic is equal to the \code{i}-the component of \code{sqrt(n)*(mpv(yuima,2)-mpv(yuima,c(1,1)))/sqrt(vartheta*mpv(yuima,r))} when \code{type="standard"}, \code{sqrt(n)*log(mpv(yuima,2)/mpv(yuima,c(1,1)))/sqrt(vartheta*mpv(yuima,r)/mpv(yuima,c(1,1))^2)} when \code{type="log"} and \code{sqrt(n)*(1-mpv(yuima,c(1,1))/mpv(yuima,2))/sqrt(vartheta*mpv(yuima,r)/mpv(yuima,c(1,1))^2)} when \code{type="ratio"}. Here, \code{n} is equal to the length of the \code{i}-th component of the \code{zoo.data} of \code{yuima} minus 1 and \code{vartheta} is \code{pi^2/4+pi-5}. When \code{adj=TRUE}, \code{mpv(yuima,r)[i]/mpv(yuima,c(1,1))^2)[i]} is replaced with 1 if it is less than 1.
+}
+\value{
+%% ~Describe the value returned
+%% If it is a LIST, use
+%% \item{comp1 }{Description of 'comp1'}
+%% \item{comp2 }{Description of 'comp2'}
+%% ...
+A list with the same length as the \code{zoo.data} of \code{yuima}. Each components of the list has class \dQuote{\code{htest}} and contains the following components:
+\item{statistic}{the value of the test statistic of the corresponding component of the \code{zoo.data} of \code{yuima}.}
+\item{p.value}{an approximate p-value for the test of the corresponding component.}
+\item{method}{the character string \dQuote{\code{Barndorff-Nielsen and Shephard jump test}}.}
+}
+\references{
+%% ~put references to the literature/web site here ~
+Barndorff-Nielsen, O. E. and Shephard, N. (2004)
+ Power and bipower variation with stochastic volatility and jumps,
+ \emph{Journal of Financial Econometrics}, \bold{2}, no. 1, 1--37.
+
+Barndorff-Nielsen, O. E. and Shephard, N. (2006)
+ Econometrics of testing for jumps in financial economics using bipower variation,
+ \emph{Journal of Financial Econometrics}, \bold{4}, no. 1, 1--30.
+
+Huang, X. and Tauchen, G. (2005)
+ The relative contribution of jumps to total price variance,
+ \emph{Journal of Financial Econometrics}, \bold{3}, no. 4, 456--499.
+}
+\author{
+%% ~~who you are~~
+The YUIMA Project Team
+}
+\note{
+%% ~~further notes~~
+Theoritically, this test may be invalid if sampling is irregular.
+}
+
+%% ~Make other sections like Warning with \section{Warning }{....} ~
+
+\seealso{
+%% ~~objects to See Also as \code{\link{help}}, ~~~
+\code{\link{mpv}}
+}
+\examples{
+
+set.seed(123)
+
+# one-dimensional case
+## Model: dXt=t*dWt+t*dzt,
+## where zt is a compound Poisson process with intensity 5 and jump sizes distribution N(0,0.1).
+
+model <- setModel(drift=0,diffusion="t",jump.coeff="t",measure.type="CP",
+ measure=list(intensity=5,df=list("dnorm(z,0,sqrt(0.1))")),
+ time.variable="t")
+
+yuima.samp <- setSampling(Terminal = 1, n = 390)
+yuima <- setYuima(model = model, sampling = yuima.samp)
+yuima <- simulate(yuima)
+plot(yuima) # The path seems to involve some jumps
+
+bns.test(yuima) # standard type
+
+bns.test(yuima,type="log") # log type
+
+bns.test(yuima,type="ratio") # ratio type
+
+# multi-dimensional case
+## Model: dXkt=t*dWk_t (k=1,2,3) (no jump case).
+
+model <- setModel(drift=c(0,0,0),diffusion=diag("t",3),time.variable="t",
+ solve.variable=c("x1","x2","x3"))
+
+yuima.samp <- setSampling(Terminal = 1, n = 390)
+yuima <- setYuima(model = model, sampling = yuima.samp)
+yuima <- simulate(yuima)
+plot(yuima)
+
+bns.test(yuima)
+
+}
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
+\keyword{ts}
+% __ONLY ONE__ keyword per line
Added: pkg/yuima/man/mpv.Rd
===================================================================
--- pkg/yuima/man/mpv.Rd (rev 0)
+++ pkg/yuima/man/mpv.Rd 2012-10-05 04:44:45 UTC (rev 206)
@@ -0,0 +1,111 @@
+\name{mpv}
+\alias{mpv,list-method}
+%- Also NEED an '\alias' for EACH other topic documented here.
+\title{
+%% ~~function to do ... ~~
+Realized Multipower variation
+}
+\description{
+%% ~~ A concise (1-5 lines) description of what the function does. ~~
+The function returns the realized MultiPower Variation (mpv), defined in Barndorff-Nielsen and Shephard (2004), for each components.
+}
+\usage{
+mpv(yuima, r = c(1, 1))
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+ \item{yuima}{
+%% ~~Describe \code{yuima} here~~
+an object of \code{\link{yuima-class}} or \code{\link{yuima.data-class}}.
+}
+ \item{r}{
+%% ~~Describe \code{r} here~~
+a vector of non-negative numbers or a list of vectors of non-negative numbers.
+}
+}
+\details{
+%% ~~ If necessary, more details than the description above ~~
+Let \code{d} be the number of the components of the \code{zoo.data} of \code{yuima}.
+
+Let \eqn{X^i_{t_0},X^i_{t_1},\dots,X^i_{t_n}} be the observation data of the \eqn{i}-th component (i.e. the \eqn{i}-th component of the \code{zoo.data} of \code{yuima}).
+
+When \eqn{r} is a \eqn{k}-dimentional vector of non-negative numbers, \code{mpv(yuima,r)} is defined as the \code{d}-dimntional vector with \code{i}-th element equal to
+\deqn{\mu_{r[1]}^{-1}\cdots\mu_{r[k]}^{-1}n^{\frac{r[1]+\cdots+r[k]}{2}-1}\sum_{j=1}^{n-k+1}|\Delta X^i_{t_{j}}|^{r[1]}|\Delta X^i_{t_{j+1}}|^{r[2]}\cdots|\Delta X^i_{t_{j+k-1}}|^{r[k]},}
+where \eqn{\mu_p} is the p-th absolute moment of the standard normal distribution and \eqn{\Delta X^i_{t_{j}}=X^i_{t_j}-X^i_{t_{j-1}}}.
+
+When \eqn{r} is a list of vectors of non-negative numbers, \code{mpv(yuima,r)} is defined as the \code{d}-dimntional vector with \code{i}-th element equal to
+\deqn{\mu_{r^i_1}^{-1}\cdots\mu_{r^i_{k_i}}^{-1}n^{\frac{r^i_1+\cdots+r^i_{k_i}}{2}-1}\sum_{j=1}^{n-k_i+1}|\Delta X^i_{t_{j}}|^{r^i_1}|\Delta X^i_{t_{j+1}}|^{r^i_2}\cdots|\Delta X^i_{t_{j+k_i-1}}|^{r^i_{k_i}},}
+where \eqn{r^i_1,\dots,r^i_{k_i}} is the \code{i}-th component of \code{r}.
+}
+\value{
+%% ~Describe the value returned
+%% If it is a LIST, use
+%% \item{comp1 }{Description of 'comp1'}
+%% \item{comp2 }{Description of 'comp2'}
+%% ...
+A numeric vector with the same length as the \code{zoo.data} of \code{yuima}
+}
+\references{
+%% ~put references to the literature/web site here ~
+Barndorff-Nielsen, O. E. and Shephard, N. (2004)
+ Power and bipower variation with stochastic volatility and jumps,
+ \emph{Journal of Financial Econometrics}, \bold{2}, no. 1, 1--37.
+
+Barndorff-Nielsen, O. E. , Graversen, S. E. , Jacod, J. , Podolskij M. and Shephard, N. (2006)
+ A central limit theorem for realised power and bipower variations of continuous semimartingales,
+ in: Kabanov, Y. , Lipster, R. , Stoyanov J. (Eds.), From Stochastic Calculus to Mathematical Finance: The Shiryaev Festschrift, Springer-Verlag, Berlin, pp. 33--68.
+}
+\author{
+%% ~~who you are~~
+The YUIMA Project Team
+}
+%\note{
+%% ~~further notes~~
+%}
+
+%% ~Make other sections like Warning with \section{Warning }{....} ~
+
+\seealso{
+%% ~~objects to See Also as \code{\link{help}}, ~~~
+\code{\link{setModel}},\code{\link{cce}}
+}
+\examples{
+
+set.seed(123)
+
+# one-dimensional case
+## Model: dXt=t*dWt+t*dzt,
+## where zt is a compound Poisson process with intensity 5 and jump sizes distribution N(0,0.1).
+
+model <- setModel(drift=0,diffusion="t",jump.coeff="t",measure.type="CP",
+ measure=list(intensity=5,df=list("dnorm(z,0,sqrt(0.1))")),
+ time.variable="t")
+
+yuima.samp <- setSampling(Terminal = 1, n = 390)
+yuima <- setYuima(model = model, sampling = yuima.samp)
+yuima <- simulate(yuima)
+plot(yuima)
+
+mpv(yuima) # true value is 1/3
+mpv(yuima,1) # true value is 1/2
+mpv(yuima,rep(2/3,3)) # true value is 1/3
+
+# multi-dimensional case
+## Model: dXkt=t*dWk_t (k=1,2,3).
+
+model <- setModel(drift=c(0,0,0),diffusion=diag("t",3),time.variable="t",
+ solve.variable=c("x1","x2","x3"))
+
+yuima.samp <- setSampling(Terminal = 1, n = 390)
+yuima <- setYuima(model = model, sampling = yuima.samp)
+yuima <- simulate(yuima)
+plot(yuima)
+
+mpv(yuima,list(c(1,1),1,rep(2/3,3))) # true varue is c(1/3,1/2,1/3)
+
+}
+
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
+\keyword{ts}
+%\keyword{ ~kwd2 }% __ONLY ONE__ keyword per line
More information about the Yuima-commits
mailing list