[Yuima-commits] r91 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jul 9 11:53:21 CEST 2010
Author: iacus
Date: 2010-07-09 11:53:21 +0200 (Fri, 09 Jul 2010)
New Revision: 91
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NEWS
pkg/yuima/R/adaBayes.R
pkg/yuima/R/qmle.R
pkg/yuima/man/quasi-likelihood.Rd
Log:
modified qmle
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2010-07-07 21:23:07 UTC (rev 90)
+++ pkg/yuima/DESCRIPTION 2010-07-09 09:53:21 UTC (rev 91)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package (unstable version)
-Version: 0.0.93
-Date: 2010-07-07
+Version: 0.0.94
+Date: 2010-07-09
Depends: methods, zoo, stats4, utils
Suggests: adapt, mvtnorm
Author: YUIMA Project Team.
Modified: pkg/yuima/NEWS
===================================================================
--- pkg/yuima/NEWS 2010-07-07 21:23:07 UTC (rev 90)
+++ pkg/yuima/NEWS 2010-07-09 09:53:21 UTC (rev 91)
@@ -1,3 +1,12 @@
+Changes in Version 0.0.94
+
+ o implemented two stage estimation in qmle
+
+ o optimized qmle for speed
+
+ o changed interface to qmle: lower and upper must
+ be a named list now
+
Changes in Version 0.0.91
o added toLatex method for yuima objects.
Modified: pkg/yuima/R/adaBayes.R
===================================================================
--- pkg/yuima/R/adaBayes.R 2010-07-07 21:23:07 UTC (rev 90)
+++ pkg/yuima/R/adaBayes.R 2010-07-09 09:53:21 UTC (rev 91)
@@ -468,7 +468,7 @@
for(j in 1:(length(slot(yuima at model@parameter,term)))){
nparam[slot(yuima at model@parameter,term)][j] <- matparam[k,][grep(slot(yuima at model@parameter,term)[j],names(unlist(param[slot(yuima at model@parameter,term)])))]
}
- qvec[k] <- exp(-negquasilik(yuima,param=nparam,print=print)+log(mdensity(nparam))+negquasilik(yuima,param=param,print=print)-log(mdensity(param)))
+ qvec[k] <- exp(quasilogl(yuima,param=nparam,print=print)+log(mdensity(nparam))-quasilogl(yuima,param=param,print=print)-log(mdensity(param)))
}
return(qvec)
}
@@ -493,7 +493,7 @@
qvec <- numeric(dim(matparam)[1])
for(k in 1:(dim(matparam)[1])){
nparam[slot(yuima at model@parameter,term)] <- matparam[k,]
- qvec[k] <- matparam[k,i]*exp(-negquasilik(yuima,param=nparam,print=print)+log(mdensity(nparam))+negquasilik(yuima,param=param,print=print)-log(mdensity(param)))
+ qvec[k] <- matparam[k,i]*exp(quasilogl(yuima,param=nparam,print=print)+log(mdensity(nparam))-quasilogl(yuima,param=param,print=print)-log(mdensity(param)))
}
return(qvec)
}
@@ -560,13 +560,13 @@
u <- runif(n.iter)
pparam <- param
- pql <- -negquasilik(yuima,param=pparam,print=print)+log(mdensity(pparam[[j]]))
+ pql <- quasilogl(yuima,param=pparam,print=print)+log(mdensity(pparam[[j]]))
nparam <- param
mhestimator <- 0
for(i in 1:n.iter){
nparam[[j]] <- pparam[[j]] + dh[,i]
- nql <- -negquasilik(yuima,param=nparam,print=print)+log(mdensity(nparam[[j]]))
+ nql <- quasilogl(yuima,param=nparam,print=print)+log(mdensity(nparam[[j]]))
alpha <- exp(nql-pql)
if(is.na(alpha)){alpha <- -Inf}
pparam[[j]] <- (alpha>u[i])*(nparam[[j]]-pparam[[j]])+pparam[[j]]
@@ -655,12 +655,12 @@
for(i in 1:n.iter){
tmp.param[[j]] <- state[,i]
- weight[i] <- exp(-negquasilik(yuima,param=liparam(tmp.param),print=print))*mdensity(liparam(tmp.param))/
+ weight[i] <- exp(quasilogl(yuima,param=liparam(tmp.param),print=print))*mdensity(liparam(tmp.param))/
dpropose(state[,i],propose.param[[1]],propose.param[[2]])
}
cstate <- param[[j]]
- cweight <- exp(-negquasilik(yuima,param=liparam(param),print=print))*mdensity(liparam(tmp.param))/
+ cweight <- exp(quasilogl(yuima,param=liparam(param),print=print))*mdensity(liparam(tmp.param))/
dpropose(param[[j]],propose.param[[1]],as.matrix(propose.param[[2]]))
u <- runif(n.iter)
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2010-07-07 21:23:07 UTC (rev 90)
+++ pkg/yuima/R/qmle.R 2010-07-09 09:53:21 UTC (rev 91)
@@ -88,7 +88,8 @@
### S.M.I. 22/06/2010
-qmle <- function(yuima, start, method="BFGS", fixed = list(), print=FALSE, ...){
+qmle <- function(yuima, start, method="BFGS", fixed = list(), print=FALSE,
+ lower, upper, ...){
call <- match.call()
@@ -96,8 +97,11 @@
yuima.stop("yuima object is missing.")
## param handling
- if( missing(start) )
- yuima.warn("Starting values for the parameters are missing. Using random initial values.")
+
+## 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
@@ -105,8 +109,12 @@
measure.par <- yuima at model@parameter at measure
common.par <- yuima at model@parameter at common
- if(length(common.par)>0)
- yuima.stop("Drift and diffusion parameters must be different.")
+ JointOptim <- FALSE
+ 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")
@@ -132,12 +140,36 @@
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(yuima:::onezoo(yuima)), env=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,]
@@ -147,20 +179,102 @@
f <- function(p) {
mycoef <- as.list(p)
- names(mycoef) <- nm
+ names(mycoef) <- nm[-idx.fixed]
mycoef[fixed.par] <- fixed
minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
}
- oout <- if(length(start)){
- optim(start, f, method = method, hessian = TRUE, ...)
+ oout <- NULL
+
+ if(length(start)){
+ if(JointOptim){ ### joint optimization
+ if(length(start)>1){ # multidimensional optim
+ oout <- optim(start, f, method = method, hessian = 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
+## DIFFUSION ESTIMATIOn first
+ old.fixed <- fixed
+ old.start <- start
+ new.start <- start[idx.diff] # considering only initial guess for diffusion
+ new.fixed <- fixed
+ 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$fixed <- NULL
+ mydots$fn <- as.name("f")
+ mydots$start <- NULL
+ mydots$par <- unlist(new.start)
+ mydots$hessian <- FALSE
+ mydots$upper <- unlist( upper[ nm[idx.diff] ])
+ mydots$lower <- unlist( lower[ nm[idx.diff] ])
+ if(length(mydots$par)>1)
+ oout <- do.call(optim, args=mydots)
+ else {
+ mydots$f <- mydots$fn
+ mydots$fn <- NULL
+ mydots$par <- NULL
+ mydots$hessian <- NULL
+ mydots$method <- NULL
+ opt1 <- do.call(optimize, args=mydots)
+ theta1 <- opt1$minimum
+ names(theta1) <- diff.par
+ oout <- list(par = theta1, value = opt1$objective)
+ }
+ theta1 <- oout$par
+
+## 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$fixed <- NULL
+ mydots$fn <- as.name("f")
+ mydots$start <- NULL
+ mydots$par <- unlist(new.start)
+ mydots$hessian <- TRUE
+ mydots$upper <- unlist( upper[ nm[idx.drift] ])
+ mydots$lower <- unlist( lower[ nm[idx.drift] ])
+
+ if(length(mydots$par)>1)
+ oout1 <- do.call(optim, args=mydots)
+ else {
+ mydots$f <- mydots$fn
+ mydots$fn <- NULL
+ mydots$par <- NULL
+ mydots$hessian <- NULL
+ mydots$method <- NULL
+ 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
+ oout1$par <- c(theta1, theta2)
+ oout <- oout1
+
+ } ### endif JointOptim
} else {
list(par = numeric(0L), value = f(start))
}
coef <- oout$par
- vcov <- if (length(coef))
- solve(oout$hessian)
+
+ vcov <- if (!is.null(oout$hessian))
+ solve(oout$hessian)
else matrix(numeric(0L), 0L, 0L)
min <- oout$value
Modified: pkg/yuima/man/quasi-likelihood.Rd
===================================================================
--- pkg/yuima/man/quasi-likelihood.Rd 2010-07-07 21:23:07 UTC (rev 90)
+++ pkg/yuima/man/quasi-likelihood.Rd 2010-07-09 09:53:21 UTC (rev 91)
@@ -17,7 +17,7 @@
ml.ql(yuima,theta2,theta1,h,theta2.lim=matrix(c(0,1),1,2),theta1.lim=matrix(c(0,1),1,2),print=FALSE,method,param,interval)
ql(yuima,theta2,theta1,h,print=FALSE,param)
rql(yuima,theta2,theta1,ptheta2,ptheta1,h,print=FALSE,param,prevparam)
-qmle(yuima, start, method="BFGS", fixed = list(), print=FALSE, ...)
+qmle(yuima, start, method="BFGS", fixed = list(), print=FALSE, lower, upper, ...)
quasilogl(yuima, param, print=FALSE)
}
\arguments{
@@ -32,6 +32,8 @@
\item{param}{\code{list} of parameters for the quasi loglikelihood.}
\item{interval}{}
\item{prevparam}{}
+ \item{lower}{a named list for specifying lower bounds of parameters}
+ \item{upper}{a named list for specifying upper bounds of parameters}
\item{start}{initial values to be passed to the optimizer.}
\item{fixed}{for conditional (quasi)maximum likelihood estimation.}
\item{...}{passed to \code{\link{optim}} method. See Examples.}
@@ -83,15 +85,14 @@
system.time(
-opt2 <- qmle(yuima, start=list(theta1=0.8, theta2=0.7), lower=c(.05,.05),
- upper=c(.5,.5), method="L-BFGS-B")
+opt2 <- qmle(yuima, start=list(theta1=0.8, theta2=0.7), lower=list(theta1=.05,theta2=.05),
+ upper=list(theta1=.5,theta2=.5), method="L-BFGS-B")
)
print("True param")
print("theta2 = .3, theta1 = .1")
print("ML estimator")
opt2 at coef
-
## another way of parameter specification
##interval <- list(theta2.lim=c(0,1), theta1.lim=c(0,1))
##system.time(
@@ -150,13 +151,25 @@
opt at coef
system.time(
-opt2 <- qmle(yuima, start=list(theta2.1=0.8, theta2.2=0.2, theta1.1=0.7, theta1.2=0.1),lower=rep(.1,4), upper=rep(1,4),method="L-BFGS-B")
+opt2 <- qmle(yuima, start=list(theta2.1=0.8, theta2.2=0.2, theta1.1=0.7, theta1.2=0.1),
+ lower=list(theta1.1=.1,theta1.2=.1,theta2.1=.1,theta2.2=.1),
+ upper=list(theta1.1=4,theta1.2=4,theta2.1=4,theta2.2=4), method="L-BFGS-B")
)
opt2 at coef
+
+## unconstrained optimization
+system.time(
+opt3 <- qmle(yuima, start=list(theta2.1=0.8, theta2.2=0.2, theta1.1=0.7, theta1.2=0.1))
+)
+opt3 at coef
+
QL
quasilogl(yuima, param=list(theta2.1=0.8, theta2.2=0.2, theta1.1=0.7, theta1.2=0.1))
+
+
+
## another way of parameter specification
#interval <- list(theta2.1.lim, theta2.2.lim, theta1.1.lim, theta1.2.lim)
#system.time(
More information about the Yuima-commits
mailing list