[Yuima-commits] r97 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Jul 11 13:09:40 CEST 2010
Author: iacus
Date: 2010-07-11 13:09:40 +0200 (Sun, 11 Jul 2010)
New Revision: 97
Added:
pkg/yuima/R/lse.R
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/R/qmle.R
pkg/yuima/man/quasi-likelihood.Rd
Log:
added lse
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2010-07-11 10:02:36 UTC (rev 96)
+++ pkg/yuima/DESCRIPTION 2010-07-11 11:09:40 UTC (rev 97)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project package (unstable version)
-Version: 0.0.97
+Version: 0.0.98
Date: 2010-07-11
Depends: methods, zoo, stats4, utils
Suggests: cubature, mvtnorm
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2010-07-11 10:02:36 UTC (rev 96)
+++ pkg/yuima/NAMESPACE 2010-07-11 11:09:40 UTC (rev 97)
@@ -86,6 +86,7 @@
export(asymptotic_term)
export(LSE)
+export(lse)
export(qmle)
export(quasilogl)
Added: pkg/yuima/R/lse.R
===================================================================
--- pkg/yuima/R/lse.R (rev 0)
+++ pkg/yuima/R/lse.R 2010-07-11 11:09:40 UTC (rev 97)
@@ -0,0 +1,96 @@
+##estimate theta2 by LSE
+##function name LSE1
+
+
+lse <- function(yuima, start, lower, upper, method="BFGS", ...){
+
+ 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
+
+ fullcoef <- c(diff.par, drift.par)
+ npar <- length(fullcoef)
+
+
+ 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)
+
+
+ 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)
+ 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)
+
+##objective function
+ f <-function(theta){
+ names(theta) <- drift.par
+ tmp <- env$deltaX - env$h * drift.term(yuima, theta, env)[-n,]
+ ret <- t(tmp) %*% tmp
+ return(sum(ret))
+ }
+
+ if(length(start)>1){ # multidimensional optim
+ oout <- optim(start, f, method = method, hessian = FALSE, lower=lower, upper=upper)
+ } else { ### one dimensional optim
+ opt1 <- optimize(f, ...) ## an interval should be provided
+ oout <- list(par = opt1$minimum, value = opt1$objective)
+ }
+
+ return(oout$par)
+}
+
+
+
+
+
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2010-07-11 10:02:36 UTC (rev 96)
+++ pkg/yuima/R/qmle.R 2010-07-11 11:09:40 UTC (rev 97)
@@ -415,7 +415,7 @@
yuima.warn(sprintf("NEG-QL: %f, %s", -QL, paste(names(param),param,sep="=",collapse=", ")))
}
if(is.infinite(QL)) return(1e10)
- return(-QL)
+ return(as.numeric(-QL))
}
Modified: pkg/yuima/man/quasi-likelihood.Rd
===================================================================
--- pkg/yuima/man/quasi-likelihood.Rd 2010-07-11 10:02:36 UTC (rev 96)
+++ pkg/yuima/man/quasi-likelihood.Rd 2010-07-11 11:09:40 UTC (rev 97)
@@ -6,13 +6,15 @@
\alias{LSE}
%\alias{ml.ql2}
\alias{rql}
+\alias{lse}
%\alias{ql,ANY-method}
%\alias{ml.ql,ANY-method}
%\alias{ml.ql2,ANY-method}
%\alias{rql,ANY-method}
-\title{Calculate quasi-likelihood and ML estimator}
+\title{Calculate quasi-likelihood and ML estimator of least squares estimator}
\description{Calculate the quasi-likelihood and estimate of the parameters of the
- stochastic differential equation by the maximum likelihood method.}
+ stochastic differential equation by the maximum likelihood method or least squares estimator
+ of the drift parameter.}
\usage{
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)
@@ -46,9 +48,15 @@
\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{lse} calculates least squares estimators of the drift parameters. This is
+ useful for initial guess of \code{qmle} estimation.
\code{quasilogl} returns the valueof the quasi loglikelihood for a given
\code{yuima} object and list of parameters \code{coef}.
+
+ The argument \code{...} can accept \code{interval} for one-dimensional
+ optimization.
}
\value{
\item{QL}{a real value.}
@@ -79,39 +87,43 @@
system.time(
opt <- ml.ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)), c(0, 1), c(0, 1))
)
-print("True param")
-print("theta2 = .3, theta1 = .1")
-print("ML estimator")
-opt at coef
+cat(sprintf("\nTrue param. theta2 = .3, theta1 = .1\n"))
+print(coef(opt))
system.time(
-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")
+opt2 <- qmle(yuima, start=list(theta1=0.8, theta2=0.7), lower=list(theta1=0,theta2=0),
+ upper=list(theta1=1,theta2=1), method="L-BFGS-B")
)
-print("True param")
-print("theta2 = .3, theta1 = .1")
-print("ML estimator")
-opt2 at coef
+cat(sprintf("\nTrue param. theta2 = .3, theta1 = .1\n"))
+print(coef(opt2))
+## initial guess for theta2 by least squares estimator
+tmp <- lse(yuima, start=list(theta2=0.7), interval=c(0,2))
+tmp
+
+system.time(
+opt3 <- qmle(yuima, start=list(theta1=0.8, theta2=tmp), lower=list(theta1=0,theta2=0),
+ upper=list(theta1=1,theta2=1), method="L-BFGS-B")
+)
+cat(sprintf("\nTrue param. theta2 = .3, theta1 = .1\n"))
+print(coef(opt3))
+
+
## perform joint estimation? Non-optimal, just for didactic purposes
system.time(
-opt3 <- 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", joint=TRUE)
+opt4 <- qmle(yuima, start=list(theta1=0.8, theta2=0.7), lower=list(theta1=0,theta2=0),
+ upper=list(theta1=1,theta2=1), method="L-BFGS-B", joint=TRUE)
)
-print("True param")
-print("theta2 = .3, theta1 = .1")
-print("ML estimator")
-opt3 at coef
+cat(sprintf("\nTrue param. theta2 = .3, theta1 = .1\n"))
+print(coef(opt4))
system.time(
-opt <- ml.ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)), c(0, 1), c(0, 1), method="Newtoon")
+opt <- ml.ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)), c(0, 1), c(0, 1), method="Newton")
)
-print("True param")
-print("theta2 = .3, theta1 = .1")
-print("ML estimator")
-opt at coef
+cat(sprintf("\nTrue param. theta2 = .3, theta1 = .1\n"))
+print(coef(opt))
##multidimension case
##dXt^e = - drift.matrix * Xt^e * dt + diff.matrix * dWt
@@ -158,12 +170,14 @@
upper=list(theta1.1=4,theta1.2=4,theta2.1=4,theta2.2=4), method="L-BFGS-B")
)
opt2 at coef
+summary(opt2)
## 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
+summary(opt3)
QL
quasilogl(yuima, param=list(theta2.1=0.8, theta2.2=0.2, theta1.1=0.7, theta1.2=0.1))
@@ -172,12 +186,6 @@
-## another way of parameter specification
-#interval <- list(theta2.1.lim, theta2.2.lim, theta1.1.lim, theta1.2.lim)
-#system.time(
-#opt <- ml.ql(yuima, h=1/((n)^(2/3)), param=param, interval=interval)
-#)
-#opt at coef
system.time(
opt <- ml.ql(yuima, theta2, theta1, h=1/((n)^(2/3)), theta2.lim, theta1.lim, method="Newton")
More information about the Yuima-commits
mailing list