[Yuima-commits] r121 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Aug 4 10:36:34 CEST 2010
Author: iacus
Date: 2010-08-04 10:36:34 +0200 (Wed, 04 Aug 2010)
New Revision: 121
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NEWS
pkg/yuima/R/CPoint.R
pkg/yuima/man/CPoint.Rd
Log:
added cpoint with symmetrized QML
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2010-08-03 18:23:47 UTC (rev 120)
+++ pkg/yuima/DESCRIPTION 2010-08-04 08:36:34 UTC (rev 121)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package (unstable version)
-Version: 0.1.11
-Date: 2010-08-03
+Version: 0.1.12
+Date: 2010-08-04
Depends: methods, zoo, stats4, utils
Suggests: cubature, mvtnorm
Author: YUIMA Project Team.
Modified: pkg/yuima/NEWS
===================================================================
--- pkg/yuima/NEWS 2010-08-03 18:23:47 UTC (rev 120)
+++ pkg/yuima/NEWS 2010-08-04 08:36:34 UTC (rev 121)
@@ -1,3 +1,7 @@
+Changes in Version 0.1.12
+
+ o added change point estimator with symmetrized QML approximation
+
Changes in Version 0.1.11
o time.variable was not copied in the env where drift and diffusion are
Modified: pkg/yuima/R/CPoint.R
===================================================================
--- pkg/yuima/R/CPoint.R 2010-08-03 18:23:47 UTC (rev 120)
+++ pkg/yuima/R/CPoint.R 2010-08-04 08:36:34 UTC (rev 121)
@@ -34,7 +34,7 @@
-CPoint <- function(yuima, param1, param2, print=FALSE, plot=FALSE){
+CPointOld <- function(yuima, param1, param2, print=FALSE, plot=FALSE){
d.size <- yuima at model@equation.number
n <- length(yuima)[1]
@@ -220,6 +220,125 @@
}
+# symmetrized version
+pminusquasiloglsym <- function(yuima, param, print=FALSE, env){
+
+ diff.par <- yuima at model@parameter at diffusion
+ fullcoef <- diff.par
+ npar <- length(fullcoef)
+
+ nm <- names(param)
+ oo <- match(nm, fullcoef)
+ if(any(is.na(oo)))
+ yuima.stop("some named arguments in 'param' are not arguments to the supplied yuima model")
+
+
+ param <- param[order(oo)]
+ nm <- names(param)
+
+ idx.diff <- match(diff.par, nm)
+
+ h <- env$h
+
+ theta1 <- unlist(param[idx.diff])
+ n.theta1 <- length(theta1)
+ n.theta <- n.theta1
+
+ d.size <- yuima at model@equation.number
+#n <- length(yuima)[1]
+ n <- dim(env$X)[1]
+
+ idx0 <- 1:round((n-1)/2)
+ vec <- matrix((env$deltaX[2*idx0+1]-2* env$deltaX[2*idx0]+env$deltaX[2*idx0-1])/sqrt(2), length(idx0), dim(env$X)[2])
+
+ K <- -0.5*d.size * log( (2*pi*h) )
+
+ print(length(idx0))
+print(dim(vec))
+
+ QL <- 0
+ pn <- numeric(length(idx0)-1)
+ diff <- diffusion.term(yuima, param, env)
+ dimB <- dim(diff[, , 1])
+
+ if(is.null(dimB)){ # one dimensional X
+ for(t in idx0[-length(idx0)]){
+ yB <- diff[, , 2*t-1]^2
+ logdet <- log(yB)
+# print(t)
+# print(vec[t, ]^2/(h*yB))
+ pn[t] <- K - 0.5*logdet-0.5*vec[t, ]^2/(h*yB)
+ QL <- QL+pn[t]
+
+ }
+ } else { # multidimensional X
+ for(t in idx0[-length(idx0)]){
+ yB <- diff[, , 2*t-1] %*% t(diff[, , 2*t-1])
+ logdet <- log(det(yB))
+ if(is.infinite(logdet) ){ # should we return 1e10?
+ pn[t] <- log(1)
+ yuima.warn("singular diffusion matrix")
+ return(1e10)
+ }else{
+ pn[t] <- K - 0.5*logdet +
+ ((-1/(2*h))*t(vec[t, ])%*%solve(yB)%*%vec[t, ])
+ QL <- QL+pn[t]
+ }
+ }
+ }
+
+ if(!is.finite(QL)){
+ yuima.warn("quasi likelihood is too small to calculate.")
+ QL <- 1e10
+ }
+
+ if(print==TRUE){
+ yuima.warn(sprintf("NEG-QL: %f, %s", -QL, paste(names(param),param,sep="=",collapse=", ")))
+ }
+
+ return(-pn)
+
+}
+
+CPoint <- function(yuima, param1, param2, print=FALSE, symmetrized=FALSE, plot=FALSE){
+ d.size <- yuima at model@equation.number
+ n <- length(yuima)[1]
+
+ d.size <- yuima at model@equation.number
+
+ env <- new.env()
+ assign("X", as.matrix(onezoo(yuima)), env=env)
+ assign("deltaX", matrix(0, n-1, d.size), env=env)
+ for(t in 1:(n-1))
+ env$deltaX[t,] <- env$X[t+1,] - env$X[t,]
+
+ assign("h", deltat(yuima at data@zoo.data[[1]]), env=env)
+ assign("time", as.numeric(index(yuima at data@zoo.data[[1]])), env=env)
+
+ QL1 <- NULL
+ QL2 <- NULL
+
+ if(!symmetrized){
+ QL1 <- pminusquasilogl(yuima=yuima, param=param1, print=print, env)
+ QL2 <- pminusquasilogl(yuima=yuima, param=param2, print=print, env)
+ } else {
+ QL1 <- pminusquasiloglsym(yuima=yuima, param=param1, print=print, env)
+ QL2 <- pminusquasiloglsym(yuima=yuima, param=param2, print=print, env)
+ }
+ nn <- length(QL1)
+ D <- sapply(2:(nn-1), function(x) sum(QL1[1:x]) + sum(QL2[-(1:x)]))
+ D <- c(D[1], D, D[length(D)])
+ if(symmetrized)
+ D <- ts(D, start=0, deltat=2*deltat(yuima at data@zoo.data[[1]]))
+ else
+ D <- ts(D, start=0, deltat=deltat(yuima at data@zoo.data[[1]]))
+ if(plot)
+ plot(D,type="l", main="change point statistics")
+# tau.hat <- index(yuima at data@zoo.data[[1]])[which.min(D)]
+ tau.hat <- index(D)[which.min(D)]
+
+ return(list(tau=tau.hat, param1=param1, param2=param2))
+}
Modified: pkg/yuima/man/CPoint.Rd
===================================================================
--- pkg/yuima/man/CPoint.Rd 2010-08-03 18:23:47 UTC (rev 120)
+++ pkg/yuima/man/CPoint.Rd 2010-08-04 08:36:34 UTC (rev 121)
@@ -5,7 +5,7 @@
\title{Adaptive LASSO estimation for stochastic differential equations}
\description{Adaptive LASSO estimation for stochastic differential equations.}
\usage{
-CPoint(yuima, param1, param2, print=FALSE, plot=FALSE)
+CPoint(yuima, param1, param2, print=FALSE, symmetrized=FALSE, plot=FALSE)
qmleL(yuima, t, ...)
qmleR(yuima, t, ...)
}
@@ -16,6 +16,7 @@
\item{plot}{plot test statistics? Default is \code{FALSE}.}
\item{print}{print some debug output. Default is \code{FALSE}.}
\item{t}{time value. See Details.}
+ \item{symmetrized}{if \code{TRUE} uses the symmetrized version of the quasi maximum-likelihood approximation.}
\item{...}{passed to \code{\link{qmle}} method. See Examples.}
}
\details{
More information about the Yuima-commits
mailing list