[Yuima-commits] r81 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jun 23 02:29:17 CEST 2010
Author: iacus
Date: 2010-06-23 02:29:17 +0200 (Wed, 23 Jun 2010)
New Revision: 81
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/R/quasi-likelihood.R
pkg/yuima/R/yuima.R
pkg/yuima/man/quasi-likelihood.Rd
Log:
added qmle and negloglik
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2010-06-20 22:51:30 UTC (rev 80)
+++ pkg/yuima/DESCRIPTION 2010-06-23 00:29:17 UTC (rev 81)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package (unstable version)
-Version: 0.0.87
-Date: 2010-06-20
+Version: 0.0.88
+Date: 2010-06-22
Depends: methods, zoo, stats4
Suggests: adapt
Author: YUIMA Project Team.
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2010-06-20 22:51:30 UTC (rev 80)
+++ pkg/yuima/NAMESPACE 2010-06-23 00:29:17 UTC (rev 81)
@@ -83,4 +83,7 @@
export(LSE)
+export(qmle)
+export(negquasilik)
+
Modified: pkg/yuima/R/quasi-likelihood.R
===================================================================
--- pkg/yuima/R/quasi-likelihood.R 2010-06-20 22:51:30 UTC (rev 80)
+++ pkg/yuima/R/quasi-likelihood.R 2010-06-23 00:29:17 UTC (rev 81)
@@ -2,6 +2,75 @@
##::extract drift term from yuima
##::para: parameter of drift term (theta2)
+
+### TO BE FIXED: all caculations should be made on a private environment to
+### avoid problems.
+### I have rewritten drift.term and diff.term instead of calc.drift and
+### calc.diffusion to make them independent of the specification of the
+### parameters. S.M.I. 22/06/2010
+
+drift.term <- function(yuima, theta){
+ r.size <- yuima at model@noise.number
+ d.size <- yuima at model@equation.number
+ modelstate <- yuima at model@state.variable
+# modelpara <- yuima at model@parameter at drift
+ DRIFT <- yuima at model@drift
+ n <- length(yuima)[1]
+ drift <- matrix(0,n,d.size)
+ X <- as.matrix(onezoo(yuima))
+# for(i in 1:length(modelpara)){
+# assign(modelpara[i],para[i])
+# }
+
+ for(i in 1:length(theta)){
+ assign(names(theta)[i],theta[[i]])
+ }
+ for(t in 1:n){
+ Xt <- X[t,]
+ for(d in 1:d.size){
+ assign(modelstate[d],Xt[d])
+ }
+ for(d in 1:d.size){
+ drift[t,d] <- eval(DRIFT[d])
+ }
+ }
+ return(drift)
+}
+
+
+
+diffusion.term <- function(yuima, theta){
+ r.size <- yuima at model@noise.number
+ d.size <- yuima at model@equation.number
+ modelstate <- yuima at model@state.variable
+# modelpara <- yuima at model@parameter at diffusion
+ DIFFUSION <- yuima at model@diffusion
+ n <- length(yuima)[1]
+ diff <- array(0, dim=c(d.size, r.size, n))
+ X <- as.matrix(onezoo(yuima))
+# for(k in 1:length(modelpara)){
+# assign(modelpara[k], para[k])
+# }
+ for(i in 1:length(theta)){
+ assign(names(theta)[i],theta[[i]])
+ }
+
+ for(r in 1:r.size){
+ for(t in 1:n){
+ Xt <- X[t, ]
+ for(d in 1:d.size){
+ assign(modelstate[d], Xt[d])
+ }
+ for(d in 1:d.size){
+ diff[d, r, t] <- eval(DIFFUSION[[d]][r])
+ }
+ }
+ }
+ return(diff)
+}
+
+
+
"calc.drift" <- function(yuima, para){
r.size <- yuima at model@noise.number
d.size <- yuima at model@equation.number
@@ -1070,3 +1139,158 @@
return(opt)
})
+
+
+
+### I have rewritten qmle as a version of ml.ql
+### This function has an interface more similar to mle.
+### ml.ql is limited in that it uses fixed names for drift and diffusion
+### parameters, while yuima model allows for any names.
+### also, I am using the same interface of optim to specify upper and lower bounds
+### S.M.I. 22/06/2010
+
+
+qmle <- function(yuima, start, method="BFGS", fixed = list(), print=FALSE, ...){
+
+ call <- match.call()
+
+ if( missing(yuima))
+ yuima.stop("yuima object is missing.")
+
+## param handling
+ if( missing(start) )
+ yuima.warn("Starting values for the parameters are missing. Using random initial values.")
+
+ 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
+
+ if(length(common.par)>0)
+ yuima.stop("Drift and diffusion parameters must be different.")
+
+ 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 <- c(diff.par, 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)
+
+ f <- function(p) {
+ mycoef <- as.list(p)
+ names(mycoef) <- nm
+ mycoef[fixed.par] <- fixed
+ negquasilik(yuima=yuima, param=mycoef, print=print)
+ }
+
+ oout <- if(length(start)){
+ optim(start, f, method = method, hessian = TRUE, ...)
+ } else {
+ list(par = numeric(0L), value = f(start))
+ }
+
+ coef <- oout$par
+ vcov <- if (length(coef))
+ solve(oout$hessian)
+ else matrix(numeric(0L), 0L, 0L)
+ min <- oout$value
+
+ mycoef <- as.list(coef)
+ names(mycoef) <- nm
+ mycoef[fixed.par] <- fixed
+
+ new("mle", call = call, coef = coef, fullcoef = unlist(mycoef),
+ vcov = vcov, min = min, details = oout, minuslogl = negquasilik,
+ method = method)
+}
+
+negquasilik <- function(yuima, param, print=FALSE){
+
+ diff.par <- yuima at model@parameter at diffusion
+ drift.par <- yuima at model@parameter at drift
+ fullcoef <- c(diff.par, 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)
+
+### FIXME
+ h <- deltat(yuima at data@zoo.data[[1]])
+
+ 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]
+ X <- as.matrix(onezoo(yuima))
+ deltaX <- matrix(0, n-1, d.size)
+ for(t in 1:(n-1))
+ deltaX[t, ] <- X[t+1, ]-X[t, ]
+
+
+ drift <- drift.term(yuima, param)
+ diff <- diffusion.term(yuima, param)
+
+ B <- calc.B(diff)
+
+ QL <- 0
+ pn <- numeric(n-1)
+ for(t in 1:(n-1)){
+ if(det(as.matrix(B[, , t]))==0){
+ pn[t] <- log(1)
+ }else{
+ pn[t] <- log( 1/((2*pi*h)^(d.size/2)*det(as.matrix(B[, , t]))^(1/2)) *
+ exp((-1/(2*h))*t(deltaX[t, ]-h*drift[t, ])%*%solve(as.matrix(B[, , t]))%*%(deltaX[t,]-h*drift[t, ])) )
+ QL <- QL+pn[t]
+ if(pn[t]==-Inf && FALSE){
+ cat("t:", t, "\n")
+ cat("B[, , t]:", B[, , t], "\n")
+ cat("det(B):",det(as.matrix(B[, , t])), "\n")
+ cat("deltaX[t, ]", deltaX[t, ], "\n")
+ cat("drift[t, ]", drift[t, ], "\n")
+ }
+ }
+ }
+ if(QL==-Inf){
+ yuima.warn("quasi likelihood is too small to calculate.")
+ }
+ if(print==TRUE){
+
+ yuima.warn(sprintf("NEG-QL: %f, %s", -QL, paste(names(param),param,sep="=",collapse=", ")))
+ }
+
+ return(-QL)
+
+}
+
+
+
Modified: pkg/yuima/R/yuima.R
===================================================================
--- pkg/yuima/R/yuima.R 2010-06-20 22:51:30 UTC (rev 80)
+++ pkg/yuima/R/yuima.R 2010-06-23 00:29:17 UTC (rev 81)
@@ -1,3 +1,6 @@
+yuima.stop <- function(x)
+stop(sprintf("\nYUIMA: %s\n", x))
+
yuima.warn <- function(x)
cat(sprintf("\nYUIMA: %s\n", x))
Modified: pkg/yuima/man/quasi-likelihood.Rd
===================================================================
--- pkg/yuima/man/quasi-likelihood.Rd 2010-06-20 22:51:30 UTC (rev 80)
+++ pkg/yuima/man/quasi-likelihood.Rd 2010-06-23 00:29:17 UTC (rev 81)
@@ -1,6 +1,8 @@
\name{ql}
\alias{ql}
\alias{ml.ql}
+\alias{qmle}
+\alias{negquasilik}
\alias{LSE}
%\alias{ml.ql2}
\alias{rql}
@@ -15,6 +17,8 @@
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, ...)
+negquasilik(yuima, param, print=FALSE)
}
\arguments{
\item{yuima}{a yuima object.}
@@ -24,15 +28,24 @@
parameters. Vector can be available only if theta is a scalar.}
\item{ptheta2,ptheta1}{}
\item{print}{you can see a progress of the estimation when print is TRUE.}
- \item{method}{}
- \item{param}{}
+ \item{method}{see Details.}
+ \item{param}{\code{list} of parameters for the negative quasi loglikelihood.}
\item{interval}{}
\item{prevparam}{}
+ \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.}
}
\details{
A function ql calculate the quasi-likelihood of a time series data X with any
parameters. A function ml.pl estimates parameters of the sdeModel by
maximizing the quasi-likelihood.
+
+ \code{qmle} behaves more likely the standard \code{mle} function in \pkg{stats4} and
+ argument \code{method} is one of the methods available in \code{\link{optim}}.
+
+ \code{negquasilik} returns the valueof the negative quasi loglikelihood for a given
+ \code{yuima} object and list of parameters \code{coef}.
}
\value{
\item{QL}{a real value.}
@@ -68,6 +81,17 @@
print("ML estimator")
opt at coef
+
+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")
+)
+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(
@@ -125,6 +149,14 @@
)
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 at coef
+QL
+-negquasilik(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