[Yuima-commits] r496 - pkg/yuima/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Oct 28 15:18:18 CEST 2016
Author: lorenzo
Date: 2016-10-28 15:18:18 +0200 (Fri, 28 Oct 2016)
New Revision: 496
Modified:
pkg/yuima/R/qmle.R
pkg/yuima/R/setCarma.R
Log:
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2016-10-28 03:45:52 UTC (rev 495)
+++ pkg/yuima/R/qmle.R 2016-10-28 13:18:18 UTC (rev 496)
@@ -10,57 +10,57 @@
### parameters. S.M.I. 22/06/2010
drift.term <- function(yuima, theta, env){
- r.size <- yuima at model@noise.number
- d.size <- yuima at model@equation.number
- modelstate <- yuima at model@state.variable
- DRIFT <- yuima at model@drift
-# n <- length(yuima)[1]
- n <- dim(env$X)[1]
+ r.size <- yuima at model@noise.number
+ d.size <- yuima at model@equation.number
+ modelstate <- yuima at model@state.variable
+ DRIFT <- yuima at model@drift
+ # n <- length(yuima)[1]
+ n <- dim(env$X)[1]
- drift <- matrix(0,n,d.size)
- tmp.env <- new.env()
- assign(yuima at model@time.variable, env$time, envir=tmp.env)
+ drift <- matrix(0,n,d.size)
+ tmp.env <- new.env()
+ assign(yuima at model@time.variable, env$time, envir=tmp.env)
- for(i in 1:length(theta)){
- assign(names(theta)[i],theta[[i]], envir=tmp.env)
- }
+ for(i in 1:length(theta)){
+ assign(names(theta)[i],theta[[i]], envir=tmp.env)
+ }
- for(d in 1:d.size){
- assign(modelstate[d], env$X[,d], envir=tmp.env)
- }
- for(d in 1:d.size){
- drift[,d] <- eval(DRIFT[d], envir=tmp.env)
- }
+ for(d in 1:d.size){
+ assign(modelstate[d], env$X[,d], envir=tmp.env)
+ }
+ for(d in 1:d.size){
+ drift[,d] <- eval(DRIFT[d], envir=tmp.env)
+ }
- return(drift)
+ return(drift)
}
diffusion.term <- function(yuima, theta, env){
- r.size <- yuima at model@noise.number
- d.size <- yuima at model@equation.number
- modelstate <- yuima at model@state.variable
- DIFFUSION <- yuima at model@diffusion
-# n <- length(yuima)[1]
- n <- dim(env$X)[1]
- tmp.env <- new.env()
- assign(yuima at model@time.variable, env$time, envir=tmp.env)
- diff <- array(0, dim=c(d.size, r.size, n))
- for(i in 1:length(theta)){
- assign(names(theta)[i],theta[[i]],envir=tmp.env)
- }
+ r.size <- yuima at model@noise.number
+ d.size <- yuima at model@equation.number
+ modelstate <- yuima at model@state.variable
+ DIFFUSION <- yuima at model@diffusion
+ # n <- length(yuima)[1]
+ n <- dim(env$X)[1]
+ tmp.env <- new.env()
+ assign(yuima at model@time.variable, env$time, envir=tmp.env)
+ diff <- array(0, dim=c(d.size, r.size, n))
+ for(i in 1:length(theta)){
+ assign(names(theta)[i],theta[[i]],envir=tmp.env)
+ }
- for(d in 1:d.size){
- assign(modelstate[d], env$X[,d], envir=tmp.env)
- }
+ for(d in 1:d.size){
+ assign(modelstate[d], env$X[,d], envir=tmp.env)
+ }
- for(r in 1:r.size){
- for(d in 1:d.size){
- diff[d, r, ] <- eval(DIFFUSION[[d]][r], envir=tmp.env)
- }
- }
- return(diff)
+ for(r in 1:r.size){
+ for(d in 1:d.size){
+ diff[d, r, ] <- eval(DIFFUSION[[d]][r], envir=tmp.env)
+ }
+ }
+ return(diff)
}
@@ -68,33 +68,33 @@
##::extract jump term from yuima
##::gamma: parameter of diffusion term (theta3)
measure.term <- function(yuima, theta, env){
- r.size <- yuima at model@noise.number
- d.size <- yuima at model@equation.number
- modelstate <- yuima at model@state.variable
- n <- dim(env$X)[1]
+ r.size <- yuima at model@noise.number
+ d.size <- yuima at model@equation.number
+ modelstate <- yuima at model@state.variable
+ n <- dim(env$X)[1]
- tmp.env <- new.env()
- assign(yuima at model@time.variable, env$time, envir =tmp.env)
- JUMP <- yuima at model@jump.coeff
- measure <- array(0, dim=c(d.size, r.size, n))
- for(i in 1:length(theta)){
- assign(names(theta)[i],theta[[i]],envir=tmp.env)
- }
+ tmp.env <- new.env()
+ assign(yuima at model@time.variable, env$time, envir =tmp.env)
+ JUMP <- yuima at model@jump.coeff
+ measure <- array(0, dim=c(d.size, r.size, n))
+ for(i in 1:length(theta)){
+ assign(names(theta)[i],theta[[i]],envir=tmp.env)
+ }
- for(d in 1:d.size){
- assign(modelstate[d], env$X[,d],envir=tmp.env)
+ for(d in 1:d.size){
+ assign(modelstate[d], env$X[,d],envir=tmp.env)
+ }
+ for(r in 1:r.size){
+ #for(d.tmp in 1:d){
+ if(d.size==1){
+ measure[1,r,] <- eval(JUMP[[r]],envir=tmp.env)
+ }else{
+ for(d in 1:d.size){
+ measure[d,r,] <- eval(JUMP[[d]][r],envir=tmp.env)
+ }
}
- for(r in 1:r.size){
- #for(d.tmp in 1:d){
- if(d.size==1){
- measure[1,r,] <- eval(JUMP[[r]],envir=tmp.env)
- }else{
- for(d in 1:d.size){
- measure[d,r,] <- eval(JUMP[[d]][r],envir=tmp.env)
- }
- }
- }
- return(measure)
+ }
+ return(measure)
}
@@ -106,23 +106,23 @@
### S.M.I. 22/06/2010
is.Poisson <- function(obj){
- if(is(obj,"yuima"))
+ if(is(obj,"yuima"))
return(is(obj at model, "yuima.poisson"))
- if(is(obj,"yuima.model"))
+ if(is(obj,"yuima.model"))
return(is(obj, "yuima.poisson"))
- return(FALSE)
+ return(FALSE)
}
is.CARMA <- function(obj){
- if(is(obj,"yuima"))
+ if(is(obj,"yuima"))
return(is(obj at model, "yuima.carma"))
- if(is(obj,"yuima.model"))
+ if(is(obj,"yuima.model"))
return(is(obj, "yuima.carma"))
- return(FALSE)
+ return(FALSE)
}
qmle <- function(yuima, start, method="BFGS", fixed = list(), print=FALSE,
- lower, upper, joint=FALSE, Est.Incr="NoIncr",aggregation=TRUE, threshold=NULL,rcpp=FALSE, ...){
+ lower, upper, joint=FALSE, Est.Incr="NoIncr",aggregation=TRUE, threshold=NULL,rcpp=FALSE, ...){
if(is(yuima at model, "yuima.carma")){
NoNeg.Noise<-FALSE
cat("\nStarting qmle for carma ... \n")
@@ -130,58 +130,58 @@
if(is.CARMA(yuima)&& length(yuima at model@info at scale.par)!=0){
method<-"L-BFGS-B"
}
- call <- match.call()
+ call <- match.call()
- if( missing(yuima))
- yuima.stop("yuima object is missing.")
- if(is.COGARCH(yuima)){
- if(missing(lower))
- lower <- list()
+ if( missing(yuima))
+ yuima.stop("yuima object is missing.")
+ if(is.COGARCH(yuima)){
+ if(missing(lower))
+ lower <- list()
- if(missing(upper))
- upper <- list()
+ if(missing(upper))
+ upper <- list()
- res <- NULL
- if("grideq" %in% names(as.list(call)[-(1:2)])){
- res <- PseudoLogLik.COGARCH(yuima, start, method=method, fixed = list(),
- lower, upper, Est.Incr, call, aggregation = aggregation, ...)
- }else{
- res <- PseudoLogLik.COGARCH(yuima, start, method=method, fixed = list(),
- lower, upper, Est.Incr, call, grideq = FALSE, aggregation = aggregation,...)
- }
+ res <- NULL
+ if("grideq" %in% names(as.list(call)[-(1:2)])){
+ res <- PseudoLogLik.COGARCH(yuima, start, method=method, fixed = list(),
+ lower, upper, Est.Incr, call, aggregation = aggregation, ...)
+ }else{
+ res <- PseudoLogLik.COGARCH(yuima, start, method=method, fixed = list(),
+ lower, upper, Est.Incr, call, grideq = FALSE, aggregation = aggregation,...)
+ }
- return(res)
- }
+ return(res)
+ }
- if(is.Ppr(yuima)){
- if(missing(lower))
- lower <- list()
+ if(is.Ppr(yuima)){
+ if(missing(lower))
+ lower <- list()
- if(missing(upper))
- upper <- list()
+ if(missing(upper))
+ upper <- list()
- # res <- NULL
- # if("grideq" %in% names(as.list(call)[-(1:2)])){
- res <- quasiLogLik.Ppr(yuimaPpr = yuima, parLambda = start, method=method, fixed = list(),
- lower, upper, call, ...)
- # }else{
- # res <- PseudoLogLik.COGARCH(yuima, start, method=method, fixed = list(),
- # lower, upper, Est.Incr, call, grideq = FALSE, aggregation = aggregation,...)
- # }
+ # res <- NULL
+ # if("grideq" %in% names(as.list(call)[-(1:2)])){
+ res <- quasiLogLik.Ppr(yuimaPpr = yuima, parLambda = start, method=method, fixed = list(),
+ lower, upper, call, ...)
+ # }else{
+ # res <- PseudoLogLik.COGARCH(yuima, start, method=method, fixed = list(),
+ # lower, upper, Est.Incr, call, grideq = FALSE, aggregation = aggregation,...)
+ # }
- return(res)
- }
+ return(res)
+ }
- orig.fixed <- fixed
- orig.fixed.par <- names(orig.fixed)
- if(is.Poisson(yuima))
- threshold <- 0
-## param handling
+ orig.fixed <- fixed
+ orig.fixed.par <- names(orig.fixed)
+ if(is.Poisson(yuima))
+ threshold <- 0
+ ## 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.")
+ ## 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.")
#14/12/2013 We modify the QMLE function when the model is a Carma(p,q).
# In this case we use a two step procedure:
@@ -189,21 +189,21 @@
# Second) Using the result in Brockwell, Davis and Yang (2007) we retrieve
# the underlying Levy. The estimated increments are used to find the L?vy parameters.
-# if(is(yuima at model, "yuima.carma")){
-# yuima.warm("two step procedure for carma(p,q)")
-# return(null)
-# }
-#
+ # if(is(yuima at model, "yuima.carma")){
+ # yuima.warm("two step procedure for carma(p,q)")
+ # return(null)
+ # }
+ #
- yuima.nobs <- as.integer(max(unlist(lapply(get.zoo.data(yuima),length))-1,na.rm=TRUE))
+ yuima.nobs <- as.integer(max(unlist(lapply(get.zoo.data(yuima),length))-1,na.rm=TRUE))
- diff.par <- yuima at model@parameter at diffusion
+ diff.par <- yuima at model@parameter at diffusion
-# 24/12
+ # 24/12
if(is.CARMA(yuima) && length(diff.par)==0
- && length(yuima at model@parameter at jump)!=0){
+ && length(yuima at model@parameter at jump)!=0){
diff.par<-yuima at model@parameter at jump
- }
+ }
if(is.CARMA(yuima) && length(yuima at model@parameter at jump)!=0){
CPlist <- c("dgamma", "dexp")
@@ -218,7 +218,7 @@
if((yuima at model@info at q+1)==(yuima at model@info at q+1)){
start[["mean.noise"]]<-1
}
- # return(NULL)
+ # return(NULL)
}
}
@@ -237,8 +237,8 @@
}
-# yuima.warn("carma(p,q): the qmle for a carma(p,q) driven by a Jump process will be implemented as soon as possible ")
-# return(NULL)
+ # yuima.warn("carma(p,q): the qmle for a carma(p,q) driven by a Jump process will be implemented as soon as possible ")
+ # return(NULL)
}
# 24/12
@@ -249,12 +249,12 @@
drift.par <- yuima at model@parameter at drift
-#01/01 we introduce the new variable in order
-# to take into account the parameters in the starting conditions
+ #01/01 we introduce the new variable in order
+ # to take into account the parameters in the starting conditions
if(is.CARMA(yuima)){
#if(length(yuima at model@info at scale.par)!=0){
- xinit.par <- yuima at model@parameter at xinit
+ xinit.par <- yuima at model@parameter at xinit
#}
}
@@ -262,14 +262,14 @@
# and CARMA, jump.par only by CARMA
jump.par <- NULL
if(is.CARMA(yuima)){
- jump.par <- yuima at model@parameter at jump
- measure.par <- yuima at model@parameter at measure
+ jump.par <- yuima at model@parameter at jump
+ measure.par <- yuima at model@parameter at measure
} else {
- if(length(yuima at model@parameter at jump)!=0){
- measure.par <- yuima at model@parameter at jump
- } else {
+ if(length(yuima at model@parameter at jump)!=0){
+ measure.par <- yuima at model@parameter at jump
+ } else {
measure.par <- yuima at model@parameter at measure
- }
+ }
}
# jump.par is used for CARMA
common.par <- yuima at model@parameter at common
@@ -279,41 +279,41 @@
if(any((match(jump.par, drift.par)))){
JointOptim <- TRUE
yuima.warn("Drift and diffusion parameters must be different. Doing
- joint estimation, asymptotic theory may not hold true.")
+ joint estimation, asymptotic theory may not hold true.")
}
+ }
+
+ 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.")
+ # 24/12
+ # if(is(yuima at model, "yuima.carma")){
+ # JointOptim <- TRUE
+ # yuima.warm("Carma(p.q): The case of common parameters in Drift and Diffusion Term will be implemented as soon as possible,")
+ # #return(NULL)
+ # }
}
- 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.")
- # 24/12
-# if(is(yuima at model, "yuima.carma")){
-# JointOptim <- TRUE
-# yuima.warm("Carma(p.q): The case of common parameters in Drift and Diffusion Term will be implemented as soon as possible,")
-# #return(NULL)
-# }
- }
+ # if(!is(yuima at model, "yuima.carma")){
+ # if(length(jump.par)+length(measure.par)>0)
+ # yuima.stop("Cannot estimate the jump models, yet")
+ # }
-# if(!is(yuima at model, "yuima.carma")){
-# 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.")
- if(!is.list(start))
- yuima.stop("Argument 'start' must be of list type.")
+ fullcoef <- NULL
- fullcoef <- NULL
+ if(length(diff.par)>0)
+ fullcoef <- diff.par
- if(length(diff.par)>0)
- fullcoef <- diff.par
+ if(length(drift.par)>0)
+ fullcoef <- c(fullcoef, drift.par)
- if(length(drift.par)>0)
- fullcoef <- c(fullcoef, drift.par)
-
if(is.CARMA(yuima) &&
- (length(yuima at model@info at loc.par)!=0)){
+ (length(yuima at model@info at loc.par)!=0)){
# 01/01 We modify the code for considering
# the loc.par in yuima.carma model
fullcoef<-c(fullcoef, yuima at model@info at loc.par)
@@ -328,13 +328,20 @@
}
# if(is.CARMA(yuima) && (length(measure.par)>0)){
- fullcoef<-c(fullcoef, measure.par)
+ fullcoef<-c(fullcoef, measure.par)
#}
- npar <- length(fullcoef)
+ if(is.CARMA(yuima)){
+ if(length(yuima at model@parameter at xinit)>1){
+ fullcoef<-unique(c(fullcoef,yuima at model@parameter at xinit))
+ }
+ }
- fixed.par <- names(fixed) # We use Fixed.par when we consider a Carma with scale parameter
+ npar <- length(fullcoef)
+
+
+ fixed.par <- names(fixed) # We use Fixed.par when we consider a Carma with scale parameter
if(is.CARMA(yuima) && (length(measure.par)>0)){
fixed.carma=NULL
if(!missing(fixed)){
@@ -374,536 +381,536 @@
for( j in c(1:length(measure.par))){
- if(is.na(match(measure.par[j],names(fixed)))){
- fixed.par <- c(fixed.par,measure.par[j])
- fixed[measure.par[j]]<-start[measure.par[j]]
+ if(is.na(match(measure.par[j],names(fixed)))){
+ fixed.par <- c(fixed.par,measure.par[j])
+ fixed[measure.par[j]]<-start[measure.par[j]]
}
}
}
- if (any(!(fixed.par %in% fullcoef)))
- yuima.stop("Some named arguments in 'fixed' are not arguments to the supplied yuima model")
+ if (any(!(fixed.par %in% fullcoef)))
+ yuima.stop("Some named arguments in 'fixed' are not arguments to the supplied yuima model")
- nm <- names(start)
+ nm <- names(start)
- oo <- match(nm, fullcoef)
+ 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)
+ 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)
- # SMI-2/9/14: idx.measure for CP
- idx.measure <- match(measure.par, nm)
- #01/01
+ idx.diff <- match(diff.par, nm)
+ idx.drift <- match(drift.par, nm)
+ # SMI-2/9/14: idx.measure for CP
+ idx.measure <- match(measure.par, nm)
+ #01/01
if(is.CARMA(yuima)){
- # if(length(yuima at model@info at scale.par)!=0){
- idx.xinit <- as.integer(na.omit(match(xinit.par,nm)))# We need to add idx if NoNeg.Noise is TRUE
+ # if(length(yuima at model@info at scale.par)!=0){
+ idx.xinit <- as.integer(na.omit(match(xinit.par,nm)))# We need to add idx if NoNeg.Noise is TRUE
#}
}
idx.fixed <- match(fixed.par, nm)
orig.idx.fixed <- idx.fixed
- 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
+ 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
+ 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
- if (is.CARMA(yuima)){
- # 24/12
+ d.size <- yuima at model@equation.number
+ if (is.CARMA(yuima)){
+ # 24/12
d.size <-1
- }
- n <- length(yuima)[1]
+ }
+ n <- length(yuima)[1]
- env <- new.env()
+ env <- new.env()
- assign("X", as.matrix(onezoo(yuima)), envir=env)
- assign("deltaX", matrix(0, n-1, d.size), envir=env)
- # SMI-2/9/14: for CP
- assign("Cn.r", numeric(n-1), envir=env)
- if(length(measure.par)==0)
- threshold <- 0 # there are no jumps, we take all observations
+ assign("X", as.matrix(onezoo(yuima)), envir=env)
+ assign("deltaX", matrix(0, n-1, d.size), envir=env)
+ # SMI-2/9/14: for CP
+ assign("Cn.r", numeric(n-1), envir=env)
+ if(length(measure.par)==0)
+ threshold <- 0 # there are no jumps, we take all observations
if (is.CARMA(yuima)){
#24/12 If we consider a carma model,
# the observations are only the first column of env$X
-# assign("X", as.matrix(onezoo(yuima)), envir=env)
-# env$X<-as.matrix(env$X[,1])
-# assign("deltaX", matrix(0, n-1, d.size)[,1], envir=env)
- env$X<-as.matrix(env$X[,1])
-# env$X<-na.omit(as.matrix(env$X[,1]))
- env$deltaX<-as.matrix(env$deltaX[,1])
+ # assign("X", as.matrix(onezoo(yuima)), envir=env)
+ # env$X<-as.matrix(env$X[,1])
+ # assign("deltaX", matrix(0, n-1, d.size)[,1], envir=env)
+ env$X<-as.matrix(env$X[,1])
+ # env$X<-na.omit(as.matrix(env$X[,1]))
+ env$deltaX<-as.matrix(env$deltaX[,1])
assign("time.obs",length(env$X),envir=env)
-# env$time.obs<-length(env$X)
+ # env$time.obs<-length(env$X)
#p <-yuima at model@info at p
assign("p", yuima at model@info at p, envir=env)
assign("q", yuima at model@info at q, envir=env)
assign("V_inf0", matrix(diag(rep(1,env$p)),env$p,env$p), envir=env)
-# env$X<-as.matrix(env$X[,1])
-# env$deltaX<-as.matrix(env$deltaX[,1])
-# assign("time.obs",length(env$X), envir=env)
-# p <-yuima at model@info at p
-# assign("V_inf0", matrix(diag(rep(1,p)),p,p), envir=env)
- }
+ # env$X<-as.matrix(env$X[,1])
+ # env$deltaX<-as.matrix(env$deltaX[,1])
+ # assign("time.obs",length(env$X), envir=env)
+ # p <-yuima at model@info at p
+ # assign("V_inf0", matrix(diag(rep(1,p)),p,p), envir=env)
+ }
assign("time", as.numeric(index(yuima at data@zoo.data[[1]])), envir=env)
- for(t in 1:(n-1)){
- env$deltaX[t,] <- env$X[t+1,] - env$X[t,]
- if(!is.CARMA(yuima))
- env$Cn.r[t] <- ((sqrt( env$deltaX[t,] %*% env$deltaX[t,])) <= threshold)
- }
+ for(t in 1:(n-1)){
+ env$deltaX[t,] <- env$X[t+1,] - env$X[t,]
+ if(!is.CARMA(yuima))
+ env$Cn.r[t] <- ((sqrt( env$deltaX[t,] %*% env$deltaX[t,])) <= threshold)
+ }
- if(length(measure.par)==0)
- env$Cn.r <- rep(1, length(env$Cn.r)) # there are no jumps, we take all observations
+ if(length(measure.par)==0)
+ env$Cn.r <- rep(1, length(env$Cn.r)) # there are no jumps, we take all observations
- assign("h", deltat(yuima at data@zoo.data[[1]]), envir=env)
+ assign("h", deltat(yuima at data@zoo.data[[1]]), envir=env)
-#SMI: 2/9/214 jump
-if(length(measure.par)>0){
+ #SMI: 2/9/214 jump
+ if(length(measure.par)>0){
args <- unlist(strsplit(suppressWarnings(sub("^.+?\\((.+)\\)", "\\1",yuima at model@measure$df$expr,perl=TRUE)), ","))
idx.intensity <- numeric(0)
for(i in 1:length(measure.par)){
- if(sum(grepl(measure.par[i],yuima at model@measure$intensity)))
+ if(sum(grepl(measure.par[i],yuima at model@measure$intensity)))
idx.intensity <- append(idx.intensity,i)
}
assign("idx.intensity", idx.intensity, envir=env)
assign("measure.var", args[1], envir=env)
-}
+ }
- f <- function(p) {
- mycoef <- as.list(p)
- if(!is.CARMA(yuima)){
- if(length(c(idx.fixed,idx.measure))>0) ## SMI 2/9/14
- names(mycoef) <- nm[-c(idx.fixed,idx.measure)] ## SMI 2/9/14
- else
- names(mycoef) <- nm
- } else {
- 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,rcpp=rcpp)
+ f <- function(p) {
+ mycoef <- as.list(p)
+ if(!is.CARMA(yuima)){
+ if(length(c(idx.fixed,idx.measure))>0) ## SMI 2/9/14
+ names(mycoef) <- nm[-c(idx.fixed,idx.measure)] ## SMI 2/9/14
+ else
+ names(mycoef) <- nm
+ } else {
+ 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,rcpp=rcpp)
+ }
-# SMI-2/9/14:
- fpsi <- function(p){
- mycoef <- as.list(p)
+ # SMI-2/9/14:
+ fpsi <- function(p){
+ mycoef <- as.list(p)
- idx.cont <- c(idx.diff,idx.drift)
- if(length(c(idx.fixed,idx.cont))>0)
- names(mycoef) <- nm[-c(idx.fixed,idx.cont)]
- else
- names(mycoef) <- nm
- mycoef[fixed.par] <- fixed
- # print(mycoef)
- #print(p)
- minusquasipsi(yuima=yuima, param=mycoef, print=print, env=env)
- }
+ idx.cont <- c(idx.diff,idx.drift)
+ if(length(c(idx.fixed,idx.cont))>0)
+ names(mycoef) <- nm[-c(idx.fixed,idx.cont)]
+ else
+ names(mycoef) <- nm
+ mycoef[fixed.par] <- fixed
+ # print(mycoef)
+ #print(p)
+ minusquasipsi(yuima=yuima, param=mycoef, print=print, env=env)
+ }
- fj <- function(p) {
- mycoef <- as.list(p)
- # names(mycoef) <- nm
- if(!is.CARMA(yuima)){
- idx.fixed <- orig.idx.fixed
- if(length(c(idx.fixed,idx.measure))>0) ## SMI 2/9/14
- names(mycoef) <- nm[-c(idx.fixed,idx.measure)] ## SMI 2/9/14
- else
- names(mycoef) <- nm
- } else {
- names(mycoef) <- nm
- mycoef[fixed.par] <- fixed
- }
- mycoef[fixed.par] <- fixed
- minusquasilogl(yuima=yuima, param=mycoef, print=print, env,rcpp=rcpp)
- }
+ fj <- function(p) {
+ mycoef <- as.list(p)
+ # names(mycoef) <- nm
+ if(!is.CARMA(yuima)){
+ idx.fixed <- orig.idx.fixed
+ if(length(c(idx.fixed,idx.measure))>0) ## SMI 2/9/14
+ names(mycoef) <- nm[-c(idx.fixed,idx.measure)] ## SMI 2/9/14
+ else
+ names(mycoef) <- nm
+ } else {
+ names(mycoef) <- nm
+ mycoef[fixed.par] <- fixed
+ }
+ mycoef[fixed.par] <- fixed
+ minusquasilogl(yuima=yuima, param=mycoef, print=print, env,rcpp=rcpp)
+ }
- oout <- NULL
- HESS <- matrix(0, length(nm), length(nm))
- colnames(HESS) <- nm
- rownames(HESS) <- nm
+ oout <- NULL
+ HESS <- matrix(0, length(nm), length(nm))
+ colnames(HESS) <- nm
+ rownames(HESS) <- nm
- HaveDriftHess <- FALSE
- HaveDiffHess <- FALSE
- HaveMeasHess <- FALSE
+ HaveDriftHess <- FALSE
+ HaveDiffHess <- FALSE
+ HaveMeasHess <- FALSE
- if(length(start)){
- if(JointOptim){ ### joint optimization
- old.fixed <- fixed
- new.start <- start
- old.start <- start
- if(!is.CARMA(yuima)){
- if(length(c(idx.fixed,idx.measure))>0)
- new.start <- start[-c(idx.fixed,idx.measure)] # considering only initial guess for
- }
+ if(length(start)){
+ if(JointOptim){ ### joint optimization
+ old.fixed <- fixed
+ new.start <- start
+ old.start <- start
+ if(!is.CARMA(yuima)){
+ if(length(c(idx.fixed,idx.measure))>0)
+ new.start <- start[-c(idx.fixed,idx.measure)] # considering only initial guess for
+ }
- if(length(new.start)>1){ #??multidimensional optim # Adjust lower for no negative Noise
- if(is.CARMA(yuima) && (NoNeg.Noise==TRUE))
- if(mean.noise %in% names(lower)){lower[mean.noise]<-10^-7}
- oout <- optim(new.start, fj, method = method, hessian = TRUE, lower=lower, upper=upper)
+ if(length(new.start)>1){ #??multidimensional optim # Adjust lower for no negative Noise
+ if(is.CARMA(yuima) && (NoNeg.Noise==TRUE))
+ if(mean.noise %in% names(lower)){lower[mean.noise]<-10^-7}
+ oout <- optim(new.start, fj, method = method, hessian = TRUE, lower=lower, upper=upper)
- if(is.CARMA(yuima)){
- HESS <- oout$hessian
- } else {
- HESS[names(new.start),names(new.start)] <- oout$hessian
- }
+ if(is.CARMA(yuima)){
+ HESS <- oout$hessian
+ } else {
+ HESS[names(new.start),names(new.start)] <- oout$hessian
+ }
- if(is.CARMA(yuima) && length(yuima at model@info at scale.par)!=0){
- b0<-paste0(yuima at model@info at ma.par,"0",collapse="")
- idx.b0<-match(b0,rownames(HESS))
- HESS<-HESS[-idx.b0,]
- HESS<-HESS[,-idx.b0]
- }
- if(is.CARMA(yuima) && length(yuima at model@parameter at measure)!=0){
- for(i in c(1:length(fixed.par))){
- indx.fixed<-match(fixed.par[i],rownames(HESS))
- HESS<-HESS[-indx.fixed,]
- HESS<-HESS[,-indx.fixed]
- }
- if(is.CARMA(yuima) && (NoNeg.Noise==TRUE)){
- idx.noise<-(match(mean.noise,rownames(HESS)))
- HESS<-HESS[-idx.noise,]
- HESS<-HESS[,-idx.noise]
- }
- }
- 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 )
- theta1 <- oout$par[diff.par]
- theta2 <- oout$par[drift.par]
+ if(is.CARMA(yuima) && length(yuima at model@info at scale.par)!=0){
+ b0<-paste0(yuima at model@info at ma.par,"0",collapse="")
+ idx.b0<-match(b0,rownames(HESS))
+ HESS<-HESS[-idx.b0,]
+ HESS<-HESS[,-idx.b0]
+ }
+ if(is.CARMA(yuima) && length(yuima at model@parameter at measure)!=0){
+ for(i in c(1:length(fixed.par))){
+ indx.fixed<-match(fixed.par[i],rownames(HESS))
+ HESS<-HESS[-indx.fixed,]
+ HESS<-HESS[,-indx.fixed]
+ }
+ if(is.CARMA(yuima) && (NoNeg.Noise==TRUE)){
+ idx.noise<-(match(mean.noise,rownames(HESS)))
+ HESS<-HESS[-idx.noise,]
+ HESS<-HESS[,-idx.noise]
+ }
+ }
+ 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 )
+ theta1 <- oout$par[diff.par]
+ theta2 <- oout$par[drift.par]
- } else { ### first diffusion, then drift
- theta1 <- NULL
+ } else { ### first diffusion, then drift
+ theta1 <- NULL
- old.fixed <- fixed
- old.start <- start
+ old.fixed <- fixed
+ old.start <- start
- if(length(idx.diff)>0){
-## DIFFUSION ESTIMATIOn first
- old.fixed <- fixed
- old.start <- start
- old.fixed.par <- fixed.par
- 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]
+ if(length(idx.diff)>0){
+ ## DIFFUSION ESTIMATIOn first
+ old.fixed <- fixed
+ old.start <- start
+ old.fixed.par <- fixed.par
+ 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$rcpp <- NULL #KK 08/07/16
- mydots$fixed <- NULL
- mydots$fn <- as.name("f")
- mydots$start <- NULL
- mydots$par <- unlist(new.start)
- mydots$hessian <- FALSE
- mydots$upper <- as.numeric(unlist( upper[ nm[idx.diff] ]))
- mydots$lower <- as.numeric(unlist( lower[ nm[idx.diff] ]))
- mydots$joint <- NULL # LM 08/03/16
- mydots$aggregation <- NULL # LM 08/03/16
- mydots$threshold <- NULL #SMI 2/9/14
+ mydots <- as.list(call)[-(1:2)]
+ mydots$print <- NULL
+ mydots$rcpp <- NULL #KK 08/07/16
+ mydots$fixed <- NULL
+ mydots$fn <- as.name("f")
+ mydots$start <- NULL
+ mydots$par <- unlist(new.start)
+ mydots$hessian <- FALSE
+ mydots$upper <- as.numeric(unlist( upper[ nm[idx.diff] ]))
+ mydots$lower <- as.numeric(unlist( lower[ nm[idx.diff] ]))
+ mydots$joint <- NULL # LM 08/03/16
+ mydots$aggregation <- NULL # LM 08/03/16
+ mydots$threshold <- NULL #SMI 2/9/14
- if((length(mydots$par)>1) | any(is.infinite(c(mydots$upper,mydots$lower)))){
- oout <- do.call(optim, args=mydots)
- } else {
- mydots$f <- mydots$fn
- mydots$fn <- NULL
- mydots$par <- NULL
- mydots$hessian <- NULL
- mydots$method <- NULL
- mydots$interval <- as.numeric(c(unlist(lower[diff.par]),unlist(upper[diff.par])))
+ if((length(mydots$par)>1) | any(is.infinite(c(mydots$upper,mydots$lower)))){
+ oout <- do.call(optim, args=mydots)
+ } else {
+ mydots$f <- mydots$fn
+ mydots$fn <- NULL
+ mydots$par <- NULL
+ mydots$hessian <- NULL
+ mydots$method <- 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
+ 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
- fixed <- old.fixed
- start <- old.start
- fixed.par <- old.fixed.par
+ fixed <- old.fixed
+ start <- old.start
+ fixed.par <- old.fixed.par
- } ## endif(length(idx.diff)>0)
+ } ## endif(length(idx.diff)>0)
- theta2 <- NULL
+ theta2 <- NULL
- if(length(idx.drift)>0){
-## DRIFT estimation with first state diffusion estimates
- fixed <- old.fixed
- start <- old.start
- old.fixed.par <- fixed.par
- 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)
+ if(length(idx.drift)>0){
+ ## DRIFT estimation with first state diffusion estimates
+ fixed <- old.fixed
+ start <- old.start
+ old.fixed.par <- fixed.par
+ 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]
+ names(new.start) <- nm[idx.drift]
- mydots <- as.list(call)[-(1:2)]
- mydots$print <- NULL
- mydots$rcpp <- NULL #KK 08/07/16
- mydots$fixed <- NULL
- mydots$fn <- as.name("f")
- mydots$threshold <- NULL #SMI 2/9/14
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 496
More information about the Yuima-commits
mailing list