[Yuima-commits] r89 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jul 7 23:21:38 CEST 2010
Author: iacus
Date: 2010-07-07 23:21:37 +0200 (Wed, 07 Jul 2010)
New Revision: 89
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/R/quasi-likelihood.R
pkg/yuima/man/quasi-likelihood.Rd
pkg/yuima/man/simulate.Rd
Log:
updated simulate examples
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2010-07-07 00:26:10 UTC (rev 88)
+++ pkg/yuima/DESCRIPTION 2010-07-07 21:21:37 UTC (rev 89)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package (unstable version)
-Version: 0.0.91
-Date: 2010-07-05
+Version: 0.0.93
+Date: 2010-07-07
Depends: methods, zoo, stats4, utils
Suggests: adapt, mvtnorm
Author: YUIMA Project Team.
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2010-07-07 00:26:10 UTC (rev 88)
+++ pkg/yuima/NAMESPACE 2010-07-07 21:21:37 UTC (rev 89)
@@ -88,7 +88,7 @@
export(LSE)
export(qmle)
-export(negquasilik)
+export(quasilogl)
S3method(toLatex, yuima)
Modified: pkg/yuima/R/quasi-likelihood.R
===================================================================
--- pkg/yuima/R/quasi-likelihood.R 2010-07-07 00:26:10 UTC (rev 88)
+++ pkg/yuima/R/quasi-likelihood.R 2010-07-07 21:21:37 UTC (rev 89)
@@ -3,74 +3,12 @@
##::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
@@ -1177,155 +1115,7 @@
-### 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/man/quasi-likelihood.Rd
===================================================================
--- pkg/yuima/man/quasi-likelihood.Rd 2010-07-07 00:26:10 UTC (rev 88)
+++ pkg/yuima/man/quasi-likelihood.Rd 2010-07-07 21:21:37 UTC (rev 89)
@@ -2,7 +2,7 @@
\alias{ql}
\alias{ml.ql}
\alias{qmle}
-\alias{negquasilik}
+\alias{quasilogl}
\alias{LSE}
%\alias{ml.ql2}
\alias{rql}
@@ -18,7 +18,7 @@
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)
+quasilogl(yuima, param, print=FALSE)
}
\arguments{
\item{yuima}{a yuima object.}
@@ -29,7 +29,7 @@
\item{ptheta2,ptheta1}{}
\item{print}{you can see a progress of the estimation when print is TRUE.}
\item{method}{see Details.}
- \item{param}{\code{list} of parameters for the negative quasi loglikelihood.}
+ \item{param}{\code{list} of parameters for the quasi loglikelihood.}
\item{interval}{}
\item{prevparam}{}
\item{start}{initial values to be passed to the optimizer.}
@@ -44,7 +44,7 @@
\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{quasilogl} returns the valueof the quasi loglikelihood for a given
\code{yuima} object and list of parameters \code{coef}.
}
\value{
@@ -154,7 +154,7 @@
)
opt2 at coef
QL
--negquasilik(yuima, param=list(theta2.1=0.8, theta2.2=0.2, theta1.1=0.7, theta1.2=0.1))
+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
Modified: pkg/yuima/man/simulate.Rd
===================================================================
--- pkg/yuima/man/simulate.Rd 2010-07-07 00:26:10 UTC (rev 88)
+++ pkg/yuima/man/simulate.Rd 2010-07-07 21:21:37 UTC (rev 89)
@@ -254,7 +254,7 @@
obj.sampling <- setSampling(Terminal=T, n=n)
obj.yuima <- setYuima(model=obj.model, sampling=obj.sampling)
X <- simulate(obj.yuima, xinit=xinit, true.parameter=list(theta=theta, sigma=sigma))
-X11()
+dev.new()
plot(X)
@@ -269,7 +269,7 @@
obj.sampling <- setSampling(Terminal=T, n=n)
obj.yuima <- setYuima(model=obj.model, sampling=obj.sampling)
X <- simulate(obj.yuima, xinit=xinit, true.parameter=list(theta=theta, sigma=sigma))
-X11()
+dev.new()
plot(X)
@@ -281,7 +281,7 @@
obj.sampling <- setSampling(Terminal=10, n=10000)
obj.yuima <- setYuima(model=obj.model, sampling=obj.sampling)
result <- simulate(obj.yuima)
-X11()
+dev.new()
plot(result)
##:: sample for multidimensional Levy process ("code" type)
More information about the Yuima-commits
mailing list