[Yuima-commits] r756 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Sep 24 11:57:00 CEST 2021
Author: kyuta
Date: 2021-09-24 11:57:00 +0200 (Fri, 24 Sep 2021)
New Revision: 756
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NEWS
pkg/yuima/R/qmle.R
pkg/yuima/man/qmle.Rd
Log:
fixed a bug regarding the fixed argument in qmle
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2021-08-17 06:12:40 UTC (rev 755)
+++ pkg/yuima/DESCRIPTION 2021-09-24 09:57:00 UTC (rev 756)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project Package for SDEs
-Version: 1.10.8
+Version: 1.10.9
Depends: R(>= 2.10.0), methods, zoo, stats4, utils, expm, cubature,
mvtnorm
Imports: Rcpp (>= 0.12.1), boot (>= 1.3-2), glassoFast, wavethresh,
Modified: pkg/yuima/NEWS
===================================================================
--- pkg/yuima/NEWS 2021-08-17 06:12:40 UTC (rev 755)
+++ pkg/yuima/NEWS 2021-09-24 09:57:00 UTC (rev 756)
@@ -75,4 +75,5 @@
2021/2/26: modified qmle.R and qmle.Rd
2021/3/13: fixed a bug in qmle when joint = TRUE and the parameter space is one-dimensional
2021/3/15: fixed a bug in qmle when yuima at model@parameter at measure is character(0) and yuima at model@measure.type is "CP"
-2021/8/16: fixed a bug in simulate when the jump intensity is not a parameter
\ No newline at end of file
+2021/8/16: fixed a bug in simulate when the jump intensity is not a parameter
+2021/9/24: fixed a bug regarding the "fixed" argument in qmle
\ No newline at end of file
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2021-08-17 06:12:40 UTC (rev 755)
+++ pkg/yuima/R/qmle.R 2021-09-24 09:57:00 UTC (rev 756)
@@ -1,3929 +1,3977 @@
-##::quasi-likelihood function
-
-##::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, 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]
-
- 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(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)
-}
-
-
-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)
- }
-
- 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)
-}
-
-
-## Koike's code
-##::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]
-
- 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(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)
-}
-
-
-### 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
-
-is.Poisson <- function(obj){
- if(is(obj,"yuima"))
- return(is(obj at model, "yuima.poisson"))
- if(is(obj,"yuima.model"))
- return(is(obj, "yuima.poisson"))
- return(FALSE)
-}
-
-is.CARMA <- function(obj){
- if(is(obj,"yuima"))
- return(is(obj at model, "yuima.carma"))
- if(is(obj,"yuima.model"))
- return(is(obj, "yuima.carma"))
- return(FALSE)
-}
-
-qmle <- function(yuima, start, method="L-BFGS-B", fixed = list(), print=FALSE,
- lower, upper, joint=FALSE, Est.Incr="NoIncr",aggregation=TRUE, threshold=NULL,rcpp=FALSE, ...){
- if(Est.Incr=="Carma.Inc"){
- Est.Incr<-"Incr"
- }
- if(Est.Incr=="Carma.Par"){
- Est.Incr<-"NoIncr"
- }
- if(Est.Incr=="Carma.IncPar"){
- Est.Incr<-"IncrPar"
- }
- if(is(yuima at model, "yuima.carma")){
- NoNeg.Noise<-FALSE
- cat("\nStarting qmle for carma ... \n")
- }
- if(is.CARMA(yuima)&& length(yuima at model@info at scale.par)!=0){
- method<-"L-BFGS-B"
- }
- call <- match.call()
-
- if( missing(yuima))
- yuima.stop("yuima object is missing.")
- if(is.COGARCH(yuima)){
- if(missing(lower))
- lower <- 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,...)
- }
-
- return(res)
- }
-
- if(is.PPR(yuima)){
- if(missing(lower))
- lower <- 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,...)
- # }
-
- return(res)
- }
-
- 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.")
-
- #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:
- # First) The Coefficient are obtained by QMLE computed using the Kalman Filter.
- # 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)
- # }
- #
-
- yuima.nobs <- as.integer(max(unlist(lapply(get.zoo.data(yuima),length))-1,na.rm=TRUE))
-
- diff.par <- yuima at model@parameter at diffusion
-
- # 24/12
- if(is.CARMA(yuima) && length(diff.par)==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")
- codelist <- c("rIG", "rgamma")
- if(yuima at model@measure.type=="CP"){
- tmp <- regexpr("\\(", yuima at model@measure$df$exp)[1]
- measurefunc <- substring(yuima at model@measure$df$exp, 1, tmp-1)
- if(!is.na(match(measurefunc,CPlist))){
- yuima.warn("carma(p,q): the qmle for a carma(p,q) driven by a Compound Poisson with no-negative random size")
- NoNeg.Noise<-TRUE
- # we need to add a new parameter for the mean of the Noise
- if((yuima at model@info at q+1)==(yuima at model@info at q+1)){
- start[["mean.noise"]]<-1
- }
- # return(NULL)
- }
-
- }
-
- if(yuima at model@measure.type=="code"){
- if(class(yuima at model@measure$df)=="yuima.law"){
- measurefunc <- "yuima.law"
- }
- else{
-
- tmp <- regexpr("\\(", yuima at model@measure$df$exp)[1]
- measurefunc <- substring(yuima at model@measure$df$exp, 1, tmp-1)
- }
- if(!is.na(match(measurefunc,codelist))){
- yuima.warn("carma(p,q): the qmle for a carma(p,q) driven by a non-Negative Levy will be implemented as soon as possible")
- NoNeg.Noise<-TRUE
- if((yuima at model@info at q+1)==(yuima at model@info at q+1)){
- start[["mean.noise"]]<-1
- }
- #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
- if(is.CARMA(yuima) && length(yuima at model@info at lin.par)>0){
- yuima.warn("carma(p,q): the case of lin.par will be implemented as soon as")
- return(NULL)
- }
-
-
- 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
-
- if(is.CARMA(yuima)){
- #if(length(yuima at model@info at scale.par)!=0){
- xinit.par <- yuima at model@parameter at xinit
- #}
- }
-
- # SMI-2/9/14: measure.par is used for Compound Poisson
- # 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
- } else {
- if(length(yuima at model@parameter at jump)!=0){
- measure.par <- unique(c(yuima at model@parameter at measure,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
-
- JointOptim <- joint
- if(is.CARMA(yuima) && length(yuima at model@parameter at jump)!=0){
- 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.")
- }
- }
-
- 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)+yuima at model@measure.type == "CP")
- # yuima.stop("Cannot estimate the jump models, yet")
- # }
-
-
- if(!is.list(start))
- yuima.stop("Argument 'start' must be of list type.")
-
- fullcoef <- NULL
-
- if(length(diff.par)>0)
- fullcoef <- diff.par
-
- if(length(drift.par)>0)
- fullcoef <- c(fullcoef, drift.par)
-
- if(is.CARMA(yuima) &&
- (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)
-
- }
-
- if(is.CARMA(yuima) && (NoNeg.Noise==TRUE)){
- if((yuima at model@info at q+1)==yuima at model@info at p){
- mean.noise<-"mean.noise"
- fullcoef<-c(fullcoef, mean.noise)
- }
- }
-
- # if(is.CARMA(yuima) && (yuima at model@measure.type == "CP")){
- fullcoef<-c(fullcoef, measure.par)
- #}
-
- if(is.CARMA(yuima)){
- if(length(yuima at model@parameter at xinit)>1){
- #fullcoef<-unique(c(fullcoef,yuima at model@parameter at xinit))
- condIniCarma<-!(yuima at model@parameter at xinit%in%fullcoef)
- if(sum(condIniCarma)>0){
- NamesInitial<- yuima at model@parameter at xinit[condIniCarma]
- start <- as.list(unlist(start)[!names(unlist(start))%in%(NamesInitial)])
- }
- }
- }
-
-
- npar <- length(fullcoef)
-
-
- fixed.par <- names(fixed) # We use Fixed.par when we consider a Carma with scale parameter
- fixed.carma=NULL
- if(is.CARMA(yuima) && (length(measure.par) > 0)){
- if(!missing(fixed)){
- if(names(fixed) %in% measure.par){
- idx.fixed.carma<-match(names(fixed),measure.par)
- idx.fixed.carma<-idx.fixed.carma[!is.na(idx.fixed.carma)]
- if(length(idx.fixed.carma)!=0){
- fixed.carma<-as.numeric(fixed[measure.par[idx.fixed.carma]])
- names(fixed.carma)<-measure.par[idx.fixed.carma]
- }
- }
- }
- upper.carma=NULL
- if(!missing(upper)){
- if(names(upper) %in% measure.par){
- idx.upper.carma<-match(names(upper),measure.par)
- idx.upper.carma<-idx.upper.carma[!is.na(idx.upper.carma)]
- if(length(idx.upper.carma)!=0){
- upper.carma<-as.numeric(upper[measure.par[idx.upper.carma]])
- names(upper.carma)<-measure.par[idx.upper.carma]
- }
- }
- }
- lower.carma=NULL
- if(!missing(lower)){
- if(names(lower) %in% measure.par){
- idx.lower.carma<-match(names(lower),measure.par)
- idx.lower.carma<-idx.lower.carma[!is.na(idx.lower.carma)]
- if(length(idx.lower.carma)!=0){
- lower.carma<-as.numeric(lower[measure.par[idx.lower.carma]])
- names(lower.carma)<-measure.par[idx.lower.carma]
- }
- }
- }
-
-
-
-
- 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 (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)
- # 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(is.null(fixed.carma)){
- idx.fixed <- match(fixed.par, nm)
- # }else{
- # dummynm <- nm[!(nm %in% fixed.par)]
- # idx.fixed <- match(fixed.par, dummynm)
- # }
- 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
-
- 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 <-1
- }
- n <- length(yuima)[1]
-
- 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(yuima at model@measure.type) == 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("time.obs",length(env$X),envir=env)
- # 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)
- }
- 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)
- }
-
- if(length(yuima at model@measure.type) == 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)
-
- #SMI: 2/9/214 jump
- if(length(yuima at model@measure.type) > 0 && yuima at model@measure.type == "CP"){
-
- # "yuima.law" LM 13/05/2018
-
- if(class(yuima at model@measure$df)=="yuima.law"){
- args <- yuima at model@parameter at measure
- }else{
- args <- unlist(strsplit(suppressWarnings(sub("^.+?\\((.+)\\)", "\\1",yuima at model@measure$df$expr,perl=TRUE)), ","))
- }
- idx.intensity <- numeric(0)
- if(length(measure.par) > 0){
- for(i in 1:length(measure.par)){
- 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)
- }
-
- # 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)
- }
-
-
- 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
-
-
- 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(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(fixed)>0)
- oout$par[fixed.par]<- unlist(fixed)[fixed.par]
-
- 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]
- # }
- # }
- if(is.CARMA(yuima)&& length(fixed)>0 && 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) && 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
- # YK Mar. 13, 2021: bug fixed
- #opt1 <- optimize(f, ...) ## an interval should be provided
- #oout <- list(par = opt1$minimum, value = opt1$objective)
- opt1 <- optimize(f, lower = lower[[names(new.start)]],
- upper = upper[[names(new.start)]], ...)
- oout <- list(par = opt1$minimum, value = opt1$objective)
- names(oout$par) <- names(new.start)
- } ### endif( length(start)>1 )
- theta1 <- oout$par[diff.par]
- theta2 <- oout$par[drift.par]
-
- } else { ### first diffusion, then drift
- theta1 <- NULL
-
- 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]
-
- 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)))){
- mydots$method<-method ##song
- oout <- do.call(optim, args=mydots)
- } else {
- mydots$f <- mydots$fn
- mydots$fn <- NULL
- mydots$par <- NULL
- mydots$hessian <- NULL
- mydots$interval <- as.numeric(c(unlist(lower[diff.par]),unlist(upper[diff.par])))
- mydots$lower <- NULL
- mydots$upper <- 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
-
- fixed <- old.fixed
- start <- old.start
- fixed.par <- old.fixed.par
-
- } ## endif(length(idx.diff)>0)
-
- 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)
-
- 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
-
- mydots$start <- NULL
- mydots$par <- unlist(new.start)
- mydots$hessian <- FALSE
- mydots$upper <- unlist( upper[ nm[idx.drift] ])
- mydots$lower <- unlist( lower[ nm[idx.drift] ])
- mydots$joint <- NULL # LM 08/03/16
- mydots$aggregation <- NULL # LM 08/03/16# LM 08/03/16
-
-
-
- if(length(mydots$par)>1 | any(is.infinite(c(mydots$upper,mydots$lower)))){
- if(is.CARMA(yuima)){
- if(NoNeg.Noise==TRUE){
- if((yuima at model@info at q+1)==yuima at model@info at p){
- mydots$lower[names(start["NoNeg.Noise"])]<-10^(-7)
- }
-
- }
- if(length(yuima at model@info at scale.par)!=0){
- name_b0<-paste0(yuima at model@info at ma.par,"0",collapse="")
- index_b0<-match(name_b0,nm)
- mydots$lower[index_b0]<-1
- mydots$upper[index_b0]<-1+10^(-7)
- }
- if (length(yuima at model@info at loc.par)!=0){
- mydots$upper <- unlist( upper[ nm ])
- mydots$lower <- unlist( lower[ nm ])
- idx.tot<-unique(c(idx.drift,idx.xinit))
- new.start <- start[idx.tot]
- names(new.start) <- nm[idx.tot]
- mydots$par <- unlist(new.start)
- }
- } # END if(is.CARMA)
-
- mydots$method <- method #song
-
- oout1 <- do.call(optim, args=mydots)
-
-
- # oout1 <- optim(mydots$par,f,method = "L-BFGS-B" , lower = mydots$lower, upper = mydots$upper)
- } else {
- mydots$f <- mydots$fn
- mydots$fn <- NULL
- mydots$par <- NULL
- mydots$hessian <- NULL
- mydots$method<-NULL
- mydots$interval <- as.numeric(c(lower[drift.par],upper[drift.par]))
-
- 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
- fixed <- old.fixed
- start <- old.start
- old.fixed.par <- fixed.par
- } ## endif(length(idx.drift)>0)
-
-
- oout1 <- list(par= c(theta1, theta2))
- if (! is.CARMA(yuima)){
- if(length(c(diff.par, diff.par))>0)
- names(oout1$par) <- c(diff.par,drift.par)
- }
-
-
- oout <- oout1
-
- } ### endif JointOptim
- } else {
- list(par = numeric(0L), value = f(start))
- }
-
-
- fMeas <- function(p) {
- mycoef <- as.list(p)
- # if(! is.CARMA(yuima)){
- # # names(mycoef) <- drift.par
- # mycoef[measure.par] <- coef[measure.par]
- #}
- minusquasipsi(yuima=yuima, param=mycoef, print=print, env=env)
- # minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
- }
-
-
- fDrift <- function(p) {
- mycoef <- as.list(p)
- if(! is.CARMA(yuima)){
- names(mycoef) <- drift.par
- mycoef[diff.par] <- coef[diff.par]
- }
- minusquasilogl(yuima=yuima, param=mycoef, print=print, env,rcpp=rcpp)
- }
-
- fDiff <- function(p) {
- mycoef <- as.list(p)
- if(! is.CARMA(yuima)){
- names(mycoef) <- diff.par
- mycoef[drift.par] <- coef[drift.par]
- }
- minusquasilogl(yuima=yuima, param=mycoef, print=print, env,rcpp=rcpp)
- }
-
- # coef <- oout$par
- #control=list()
- #par <- coef
-
- #names(par) <- unique(c(diff.par, drift.par))
- # nm <- unique(c(diff.par, drift.par))
-
- # START: ESTIMATION OF CP part
- theta3 <- NULL
-
- if(length(idx.measure)>0 & !is.CARMA(yuima)){
- idx.cont <- c(idx.drift,idx.diff)
-
- fixed <- old.fixed
- start <- old.start
- old.fixed.par <- fixed.par
- new.fixed <- fixed
-
- new.start <- start[idx.measure] # considering only initial guess for measure
- new.fixed <- fixed
-
- new.fixed[names(theta1)] <- theta1
- new.fixed[names(theta2)] <- theta2
-
- fixed <- new.fixed
- fixed.par <- names(fixed)
- idx.fixed <- match(fixed.par, nm)
- # names(new.start) <- nm[idx.drift]
- names(new.start) <- nm[idx.measure]
-
- mydots <- as.list(call)[-(1:2)]
- # mydots$print <- NULL
- mydots$threshold <- NULL
- mydots$fixed <- NULL
- mydots$fn <- as.name("fpsi")
- mydots$start <- NULL
- mydots$threshold <- NULL #SMI 2/9/14
-
- mydots$par <- unlist(new.start)
- mydots$hessian <- TRUE
- mydots$joint <- NULL
- mydots$upper <- unlist( upper[ nm[idx.measure] ])
- mydots$lower <- unlist( lower[ nm[idx.measure] ])
- mydots$method <- method
-
- oout3 <- do.call(optim, args=mydots)
-
- theta3 <- oout3$par
- #print(theta3)
- HESS[measure.par,measure.par] <- oout3$hessian
- HaveMeasHess <- TRUE
-
- fixed <- old.fixed
- start <- old.start
- fixed.par <- old.fixed.par
- }
- # END: ESTIMATION OF CP part
-
-
-
- if(!is.CARMA(yuima)){
-
- oout4 <- list(par= c(theta1, theta2, theta3))
- names(oout4$par) <- c(diff.par,drift.par,measure.par)
- oout <- oout4
- }
-
- coef <- oout$par
-
-
- control=list()
- par <- coef
- if(!is.CARMA(yuima)){
-
- names(par) <- unique(c(diff.par, drift.par,measure.par))
- nm <- unique(c(diff.par, drift.par,measure.par))
- } else {
- names(par) <- unique(c(diff.par, drift.par))
- nm <- unique(c(diff.par, drift.par))
- }
- #return(oout)
-
-
- if(is.CARMA(yuima) && length(yuima at model@parameter at measure)!=0){
- nm <-c(nm,measure.par)
- if((NoNeg.Noise==TRUE)){nm <-c(nm,mean.noise)}
-
- nm<-unique(nm)
- }
- if(is.CARMA(yuima) && (length(yuima at model@info at loc.par)!=0)){
- nm <-unique(c(nm,yuima at model@info at loc.par))
- }
-
-
- conDrift <- list(trace = 5, fnscale = 1,
- parscale = rep.int(5, length(drift.par)),
- ndeps = rep.int(0.001, length(drift.par)), maxit = 100L,
- abstol = -Inf, reltol = sqrt(.Machine$double.eps), alpha = 1,
- beta = 0.5, gamma = 2, REPORT = 10, type = 1, lmm = 5,
- factr = 1e+07, pgtol = 0, tmax = 10, temp = 10)
- conDiff <- list(trace = 5, fnscale = 1,
- parscale = rep.int(5, length(diff.par)),
- ndeps = rep.int(0.001, length(diff.par)), maxit = 100L,
- abstol = -Inf, reltol = sqrt(.Machine$double.eps), alpha = 1,
- beta = 0.5, gamma = 2, REPORT = 10, type = 1, lmm = 5,
- factr = 1e+07, pgtol = 0, tmax = 10, temp = 10)
- conMeas <- list(trace = 5, fnscale = 1,
- parscale = rep.int(5, length(measure.par)),
- ndeps = rep.int(0.001, length(measure.par)), maxit = 100L,
- abstol = -Inf, reltol = sqrt(.Machine$double.eps), alpha = 1,
- beta = 0.5, gamma = 2, REPORT = 10, type = 1, lmm = 5,
- factr = 1e+07, pgtol = 0, tmax = 10, temp = 10)
- if(is.CARMA(yuima) && length(yuima at model@info at loc.par)!=0 ){
- conDrift <- list(trace = 5, fnscale = 1,
- parscale = rep.int(5, length(c(drift.par,yuima at model@info at loc.par))),
- ndeps = rep.int(0.001, length(c(drift.par,yuima at model@info at loc.par))),
- maxit = 100L,
- abstol = -Inf, reltol = sqrt(.Machine$double.eps), alpha = 1,
- beta = 0.5, gamma = 2, REPORT = 10, type = 1, lmm = 5,
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 756
More information about the Yuima-commits
mailing list