[Yuima-commits] r171 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Aug 8 14:01:16 CEST 2011
Author: iacus
Date: 2011-08-08 14:01:15 +0200 (Mon, 08 Aug 2011)
New Revision: 171
Added:
pkg/yuima/R/phi.test.R
pkg/yuima/man/phi.test.Rd
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/R/qmle.R
Log:
added phi div test stat
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2011-08-04 01:40:15 UTC (rev 170)
+++ pkg/yuima/DESCRIPTION 2011-08-08 12:01:15 UTC (rev 171)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package (unstable version)
-Version: 0.1.1931
-Date: 2011-08-04
+Version: 0.1.1932
+Date: 2011-08-08
Depends: methods, zoo, stats4, utils
Suggests: cubature, mvtnorm
Author: YUIMA Project Team.
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2011-08-04 01:40:15 UTC (rev 170)
+++ pkg/yuima/NAMESPACE 2011-08-08 12:01:15 UTC (rev 171)
@@ -5,12 +5,14 @@
importFrom("zoo")
importFrom(stats, confint)
+
importFrom(utils, toLatex)
+
exportClasses("yuima",
"yuima.data",
"yuima.sampling",
@@ -94,6 +96,7 @@
export(qmle)
export(quasilogl)
+export(phi.test)
export(lasso)
export(CPoint)
export(qmleR)
@@ -101,6 +104,7 @@
export(cbind.yuima)
+S3method(print, phitest)
S3method(toLatex, yuima)
S3method(toLatex, yuima.model)
Added: pkg/yuima/R/phi.test.R
===================================================================
--- pkg/yuima/R/phi.test.R (rev 0)
+++ pkg/yuima/R/phi.test.R 2011-08-08 12:01:15 UTC (rev 171)
@@ -0,0 +1,49 @@
+
+phi.test <- function(yuima, H0, H1, phi=log, print=FALSE,...){
+
+ d.size <- yuima at model@equation.number
+ n <- length(yuima)[1]
+
+ env <- new.env()
+ assign("X", as.matrix(yuima:::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)
+ est <- FALSE
+ if(missing(H1)){
+ cat("\nestimating parameters via QMLE...\n")
+ H1 <- coef(qmle(yuima, ...))
+ est <- TRUE
+ }
+ g0 <- quasiloglvec(yuima=yuima, param=H0, print=print, env)
+ g1 <- quasiloglvec(yuima=yuima, param=H1, print=print, env)
+ div <- mean(phi(g1-g0), na.rm=TRUE)
+ stat <- 2*sum(phi(g1-g0), na.rm=TRUE)
+ df <- length(H0)
+ val <- list(div=div, stat=stat, H0=H0, H1=H1, phi=deparse(substitute(phi)), pvalue=1-pchisq(stat, df=df), df=df,est=est)
+ attr(val, "class") <- "phitest"
+ return( val )
+}
+
+
+
+print.phitest <- function(x,...){
+
+ symnum(x$pvalue, corr = FALSE, na = FALSE,
+ cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
+ symbols = c("***", "**", "*", ".", " "))->Signif
+ cat(sprintf("Phi-Divergence test statistic based on phi function '%s'\n",x$phi) )
+ nm <- names(x$H0)
+ cat("H0: ")
+ cat(sprintf("%s = %s", nm, format(x$H0,digits=3,nsmall=3)))
+ cat("\nversus\n")
+ cat("H1: ")
+ cat(sprintf("%s = %s", nm, format(x$H1[nm],digits=3,nsmall=3)))
+ if(x$est)
+ cat("\nH1 parameters estimated using QMLE")
+ cat(sprintf("\n\nTest statistic = %s, df = %d, p-value = %s %s\n", format(x$stat,digits=2,nsmall=3), x$df, format(x$pvalue), Signif))
+ cat("---\nSignif. codes: ", attr(Signif, "legend"), "\n")
+}
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2011-08-04 01:40:15 UTC (rev 170)
+++ pkg/yuima/R/qmle.R 2011-08-08 12:01:15 UTC (rev 171)
@@ -510,3 +510,422 @@
}
+
+
+
+
+# returns the vector of log-transitions instead of the final quasilog
+quasiloglvec <- function(yuima, param, print=FALSE, env){
+
+ diff.par <- yuima at model@parameter at diffusion
+ drift.par <- yuima at model@parameter at drift
+
+ fullcoef <- NULL
+
+ if(length(diff.par)>0)
+ fullcoef <- diff.par
+
+ if(length(drift.par)>0)
+ fullcoef <- c(fullcoef, drift.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)
+ idx.drift <- match(drift.par, nm)
+
+ h <- env$h
+
+ theta1 <- unlist(param[idx.diff])
+ theta2 <- unlist(param[idx.drift])
+ n.theta1 <- length(theta1)
+ n.theta2 <- length(theta2)
+ n.theta <- n.theta1+n.theta2
+
+ d.size <- yuima at model@equation.number
+ n <- length(yuima)[1]
+
+
+ drift <- drift.term(yuima, param, env)
+ diff <- diffusion.term(yuima, param, env)
+
+ QL <- numeric(n-1) ## here is the difference
+
+ pn <- 0
+
+
+ vec <- env$deltaX-h*drift[-n,]
+
+ K <- -0.5*d.size * log( (2*pi*h) )
+
+ dimB <- dim(diff[, , 1])
+
+ if(is.null(dimB)){ # one dimensional X
+ for(t in 1:(n-1)){
+ yB <- diff[, , t]^2
+ logdet <- log(yB)
+ pn <- K - 0.5*logdet-0.5*vec[t, ]^2/(h*yB)
+ QL[t] <- pn
+
+ }
+ } else { # multidimensional X
+ for(t in 1:(n-1)){
+ yB <- diff[, , t] %*% t(diff[, , t])
+ logdet <- log(det(yB))
+ if(is.infinite(logdet) ){ # should we return 1e10?
+ pn <- log(1)
+ yuima.warn("singular diffusion matrix")
+ return(1e10)
+ }else{
+ pn <- K - 0.5*logdet +
+ ((-1/(2*h))*t(vec[t, ])%*%solve(yB)%*%vec[t, ])
+ QL[t] <- pn
+ }
+ }
+ }
+ return(QL)
+}
+
+
+
+
+## test function using nlmnb instead of optim. Not mush difference
+
+qmle2 <- function(yuima, start, method="BFGS", fixed = list(), print=FALSE,
+lower, upper, joint=FALSE, ...){
+
+ call <- match.call()
+
+ if( missing(yuima))
+ yuima.stop("yuima object is missing.")
+
+ ## param handling
+
+ ## FIXME: maybe we should choose initial values at random within lower/upper
+ ## at present, qmle stops
+ if( missing(start) )
+ yuima.stop("Starting values for the parameters are missing.")
+
+ diff.par <- yuima at model@parameter at diffusion
+ drift.par <- yuima at model@parameter at drift
+ jump.par <- yuima at model@parameter at jump
+ measure.par <- yuima at model@parameter at measure
+ common.par <- yuima at model@parameter at common
+
+ JointOptim <- joint
+ if(length(common.par)>0){
+ JointOptim <- TRUE
+ yuima.warn("Drift and diffusion parameters must be different. Doing
+ joint estimation, asymptotic theory may not hold true.")
+ }
+
+
+ if(length(jump.par)+length(measure.par)>0)
+ yuima.stop("Cannot estimate the jump models, yet")
+
+ if(!is.list(start))
+ yuima.stop("Argument 'start' must be of list type.")
+
+ fullcoef <- NULL
+
+ if(length(diff.par)>0)
+ fullcoef <- diff.par
+
+ if(length(drift.par)>0)
+ fullcoef <- c(fullcoef, drift.par)
+
+ npar <- length(fullcoef)
+
+ fixed.par <- names(fixed)
+
+ if (any(!(fixed.par %in% fullcoef)))
+ yuima.stop("Some named arguments in 'fixed' are not arguments to the supplied yuima model")
+
+ nm <- names(start)
+ oo <- match(nm, fullcoef)
+ if(any(is.na(oo)))
+ yuima.stop("some named arguments in 'start' are not arguments to the supplied yuima model")
+ start <- start[order(oo)]
+ nm <- names(start)
+
+ idx.diff <- match(diff.par, nm)
+ idx.drift <- match(drift.par, nm)
+ idx.fixed <- match(fixed.par, nm)
+
+ tmplower <- as.list( rep( -Inf, length(nm)))
+ names(tmplower) <- nm
+ if(!missing(lower)){
+ idx <- match(names(lower), names(tmplower))
+ if(any(is.na(idx)))
+ yuima.stop("names in 'lower' do not match names fo parameters")
+ tmplower[ idx ] <- lower
+ }
+ lower <- tmplower
+
+ tmpupper <- as.list( rep( Inf, length(nm)))
+ names(tmpupper) <- nm
+ if(!missing(upper)){
+ idx <- match(names(upper), names(tmpupper))
+ if(any(is.na(idx)))
+ yuima.stop("names in 'lower' do not match names fo parameters")
+ tmpupper[ idx ] <- upper
+ }
+ upper <- tmpupper
+
+
+
+
+ d.size <- yuima at model@equation.number
+ n <- length(yuima)[1]
+
+ env <- new.env()
+ assign("X", as.matrix(onezoo(yuima)), env=env)
+ assign("deltaX", matrix(0, n-1, d.size), env=env)
+ assign("time", as.numeric(index(yuima at data@zoo.data[[1]])), 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)
+
+ f <- function(p) {
+ mycoef <- as.list(p)
+ # print(nm[-idx.fixed])
+ # print(nm)
+ if(length(idx.fixed)>0)
+ names(mycoef) <- nm[-idx.fixed]
+ else
+ names(mycoef) <- nm
+ mycoef[fixed.par] <- fixed
+ minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
+ }
+
+ fj <- function(p) {
+ mycoef <- as.list(p)
+ names(mycoef) <- nm
+ mycoef[fixed.par] <- fixed
+ minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
+ }
+
+ oout <- NULL
+ HESS <- matrix(0, length(nm), length(nm))
+ colnames(HESS) <- nm
+ rownames(HESS) <- nm
+ HaveDriftHess <- FALSE
+ HaveDiffHess <- FALSE
+ if(length(start)){
+ if(JointOptim){ ### joint optimization
+ if(length(start)>1){ # multidimensional optim
+ oout <- nlminb(start, fj, lower=lower, upper=upper)
+ print(oout)
+ HESS <- oout$hessian
+ HaveDriftHess <- TRUE
+ HaveDiffHess <- TRUE
+ } else { ### one dimensional optim
+ opt1 <- optimize(f, ...) ## an interval should be provided
+ oout <- list(par = opt1$minimum, value = opt1$objective)
+ } ### endif( length(start)>1 )
+ } else { ### first diffusion, then drift
+ theta1 <- NULL
+
+ old.fixed <- fixed
+ old.start <- start
+
+ if(length(idx.diff)>0){
+ ## DIFFUSION ESTIMATIOn first
+ old.fixed <- fixed
+ old.start <- start
+ new.start <- start[idx.diff] # considering only initial guess for diffusion
+ new.fixed <- fixed
+ if(length(idx.drift)>0)
+ new.fixed[nm[idx.drift]] <- start[idx.drift]
+ fixed <- new.fixed
+ fixed.par <- names(fixed)
+ idx.fixed <- match(fixed.par, nm)
+ names(new.start) <- nm[idx.diff]
+
+ mydots <- as.list(call)[-(1:2)]
+ mydots$print <- NULL
+ mydots$fixed <- NULL
+ mydots$fn <- as.name("f")
+ mydots$objective <- mydots$fn
+ mydots$start <- NULL
+ mydots$par <- unlist(new.start)
+ mydots$start <- unlist(new.start)
+ mydots$hessian <- NULL
+ mydots$upper <- unlist( upper[ nm[idx.diff] ])
+ mydots$lower <- unlist( lower[ nm[idx.diff] ])
+ if(length(mydots$par)>1){
+ mydots$fn <- NULL
+ mydots$par <- NULL
+
+ oout <- do.call(nlminb, args=mydots)
+ # oout <- do.call(optim, args=mydots)
+ } else {
+ mydots$f <- mydots$fn
+ mydots$fn <- NULL
+ mydots$par <- NULL
+ mydots$hessian <- NULL
+ mydots$method <- NULL
+ mydots$start <- NULL
+ mydots$objective <- NULL
+ mydots$interval <- as.numeric(c(unlist(lower[diff.par]),unlist(upper[diff.par])))
+ mydots$lower <- NULL
+ mydots$upper <- NULL
+ opt1 <- do.call(optimize, args=mydots)
+ theta1 <- opt1$minimum
+ names(theta1) <- diff.par
+ oout <- list(par = theta1, value = opt1$objective)
+ }
+ theta1 <- oout$par
+ } ## endif(length(idx.diff)>0)
+
+ theta2 <- NULL
+
+ if(length(idx.drift)>0){
+ ## DRIFT estimation with first state diffusion estimates
+ fixed <- old.fixed
+ start <- old.start
+ new.start <- start[idx.drift] # considering only initial guess for drift
+ new.fixed <- fixed
+ new.fixed[names(theta1)] <- theta1
+ fixed <- new.fixed
+ fixed.par <- names(fixed)
+ idx.fixed <- match(fixed.par, nm)
+ names(new.start) <- nm[idx.drift]
+
+ mydots <- as.list(call)[-(1:2)]
+ mydots$print <- NULL
+ mydots$fixed <- NULL
+ mydots$fn <- as.name("f")
+ mydots$objective <- mydots$fn
+ mydots$start <- NULL
+ mydots$par <- unlist(new.start)
+ mydots$start <- unlist(new.start)
+ mydots$hessian <- NULL
+ mydots$upper <- unlist( upper[ nm[idx.drift] ])
+ mydots$lower <- unlist( lower[ nm[idx.drift] ])
+
+ if(length(mydots$par)>1){
+ mydots$par <- NULL
+
+
+mydots$fn <- NULL
+ oout1 <- do.call(nlminb, args=mydots)
+ # oout1 <- do.call(optim, args=mydots)
+ } else {
+ mydots$f <- mydots$fn
+ mydots$fn <- NULL
+ mydots$par <- NULL
+ mydots$start <- NULL
+ mydots$objective <- NULL
+ mydots$hessian <- NULL
+ mydots$method <- NULL
+ mydots$interval <- as.numeric(c(lower[drift.par],upper[drift.par]))
+ opt1 <- do.call(optimize, args=mydots)
+ theta2 <- opt1$minimum
+ names(theta2) <- drift.par
+ oout1 <- list(par = theta2, value = as.numeric(opt1$objective))
+ }
+ theta2 <- oout1$par
+ } ## endif(length(idx.drift)>0)
+ oout1 <- list(par= c(theta1, theta2))
+ names(oout1$par) <- c(diff.par,drift.par)
+ oout <- oout1
+
+ } ### endif JointOptim
+ } else {
+ list(par = numeric(0L), value = f(start))
+ }
+
+
+ fDrift <- function(p) {
+ mycoef <- as.list(p)
+ names(mycoef) <- drift.par
+ mycoef[diff.par] <- coef[diff.par]
+ minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
+ }
+
+ fDiff <- function(p) {
+ mycoef <- as.list(p)
+ names(mycoef) <- diff.par
+ mycoef[drift.par] <- coef[drift.par]
+ minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
+ }
+
+ coef <- oout$par
+ control=list()
+ par <- coef
+ names(par) <- c(diff.par, drift.par)
+ nm <- c(diff.par, drift.par)
+
+ # print(par)
+ # print(coef)
+ conDrift <- list(trace = 5, fnscale = 1,
+ parscale = rep.int(5, length(drift.par)),
+ ndeps = rep.int(0.001, length(drift.par)), maxit = 100L,
+ abstol = -Inf, reltol = sqrt(.Machine$double.eps), alpha = 1,
+ beta = 0.5, gamma = 2, REPORT = 10, type = 1, lmm = 5,
+ factr = 1e+07, pgtol = 0, tmax = 10, temp = 10)
+ conDiff <- list(trace = 5, fnscale = 1,
+ parscale = rep.int(5, length(diff.par)),
+ ndeps = rep.int(0.001, length(diff.par)), maxit = 100L,
+ abstol = -Inf, reltol = sqrt(.Machine$double.eps), alpha = 1,
+ beta = 0.5, gamma = 2, REPORT = 10, type = 1, lmm = 5,
+ factr = 1e+07, pgtol = 0, tmax = 10, temp = 10)
+
+ # nmsC <- names(con)
+ # if (method == "Nelder-Mead")
+ # con$maxit <- 500
+ # if (method == "SANN") {
+ # con$maxit <- 10000
+ # con$REPORT <- 100
+ # }
+ # con[(namc <- names(control))] <- control
+ # if (length(noNms <- namc[!namc %in% nmsC]))
+ # warning("unknown names in control: ", paste(noNms, collapse = ", "))
+ # if (con$trace < 0)
+ # warning("read the documentation for 'trace' more carefully")
+ # else if (method == "SANN" && con$trace && as.integer(con$REPORT) ==
+ # 0)
+ # stop("'trace != 0' needs 'REPORT >= 1'")
+ # if (method == "L-BFGS-B" && any(!is.na(match(c("reltol",
+ # "abstol"), namc))))
+ # warning("method L-BFGS-B uses 'factr' (and 'pgtol') instead of 'reltol' and 'abstol'")
+ # npar <- length(par)
+ # if (npar == 1 && method == "Nelder-Mead")
+ # warning("one-diml optimization by Nelder-Mead is unreliable: use optimize")
+ #
+ if(!HaveDriftHess & (length(drift.par)>0)){
+ hess2 <- .Internal(optimhess(coef[drift.par], fDrift, NULL, conDrift))
+ HESS[drift.par,drift.par] <- hess2
+ }
+
+ if(!HaveDiffHess & (length(diff.par)>0)){
+ hess1 <- .Internal(optimhess(coef[diff.par], fDiff, NULL, conDiff))
+ HESS[diff.par,diff.par] <- hess1
+ }
+
+ oout$hessian <- HESS
+
+ vcov <- if (length(coef))
+ solve(oout$hessian)
+ else matrix(numeric(0L), 0L, 0L)
+
+ mycoef <- as.list(coef)
+ names(mycoef) <- nm
+ mycoef[fixed.par] <- fixed
+
+ min <- minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
+
+ new("mle", call = call, coef = coef, fullcoef = unlist(mycoef),
+ vcov = vcov, min = min, details = oout, minuslogl = minusquasilogl,
+ method = method)
+}
+
+
Added: pkg/yuima/man/phi.test.Rd
===================================================================
--- pkg/yuima/man/phi.test.Rd (rev 0)
+++ pkg/yuima/man/phi.test.Rd 2011-08-08 12:01:15 UTC (rev 171)
@@ -0,0 +1,46 @@
+\name{phi.test}
+\alias{phi.test}
+\title{Phi-divergence test statistic for stochastic differential equations}
+\description{Phi-divergence test statistic for stochastic differential equations.}
+\usage{
+phi.test(yuima, H0, H1, phi=log, print=FALSE,...)
+}
+\arguments{
+ \item{yuima}{a yuima object.}
+ \item{H0}{a named list of parameter under H0.}
+ \item{H1}{a named list of parameter under H1.}
+ \item{phi}{the phi function to be used in the test.}
+ \item{print}{you can see a progress of the estimation when print is \code{TRUE}.}
+ \item{...}{passed to \code{\link{qmle}} function.}
+}
+\details{
+
+ \code{phi.test} executes a Phi-divergence test. If \code{H1} is not specified
+ this hypothesis is filled with the QMLE estimates.
+
+}
+\value{
+ \item{ans}{an obkect of class \code{phitest}.}
+}
+\author{The YUIMA Project Team}
+\examples{
+model<- setModel(drift="t1*(t2-x)",diffusion="t3*x")
+T<-10
+n<-1000
+sampling <- setSampling(Terminal=T,n=n)
+yuima<-setYuima(model=model, sampling=sampling)
+
+h0 <- list(t1=0.8, t2=0.2, t3=0.5)
+X <- simulate(yuima,xinit=10,true=h0)
+h1 <- list(t1=0.3, t2=0.2, t3=0.1)
+
+phi1 <- function(x) 1-x+x*log(x)
+
+phi.test(X, H0=h0, H1=h1,phi=phi1)
+phi.test(X, H0=h0, phi=phi1, start=h0, lower=list(t1=0.1, t2=0.1, t3=0.1), upper=list(t1=2,t2=2,t3=2),method="L-BFGS-B")
+phi.test(X, H0=h1, phi=phi1, start=h0, lower=list(t1=0.1, t2=0.1, t3=0.1), upper=list(t1=2,t2=2,t3=2),method="L-BFGS-B")
+}
+
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
+\keyword{ts}
More information about the Yuima-commits
mailing list