[Yuima-commits] r745 - pkg/yuima/R

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Fri Feb 26 11:35:47 CET 2021


Author: eguchi
Date: 2021-02-26 11:35:47 +0100 (Fri, 26 Feb 2021)
New Revision: 745

Modified:
   pkg/yuima/R/qmle.R
Log:
modified

Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R	2021-02-05 04:55:07 UTC (rev 744)
+++ pkg/yuima/R/qmle.R	2021-02-26 10:35:47 UTC (rev 745)
@@ -1,3923 +1,3922 @@
-##::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="BFGS", 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)+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 <- 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) && (length(measure.par)>0)){
-  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(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("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(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)
-
-  #SMI: 2/9/214 jump
-  if(length(measure.par)>0){
-
-  #  "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)
-    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
-        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
-
-      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)))){
-          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
-
-        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)
-
-
-
-          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,
-                     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)
-  }
-
-
-
-  if(!HaveDriftHess & (length(drift.par)>0)){
-    #hess2 <- .Internal(optimhess(coef[drift.par], fDrift, NULL, conDrift))
-    if(!is.CARMA(yuima)){
-      hess2 <- optimHess(coef[drift.par], fDrift, NULL, control=conDrift)
-      HESS[drift.par,drift.par] <- hess2
-    } else{
-      names(coef) <- c(drift.par,yuima at model@info at loc.par)
-      hess2 <- optimHess(coef, fDrift, NULL, control=conDrift)
-      HESS <- hess2
-    }
-    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(!HaveDiffHess  & (length(diff.par)>0)){
-    hess1 <- optimHess(coef[diff.par], fDiff, NULL, control=conDiff)
-    HESS[diff.par,diff.par] <- hess1
-  }
-
-  oout$hessian <- HESS
-
-
-  if(!HaveMeasHess & (length(measure.par)>0) & !is.CARMA(yuima)){
-    hess1 <- optimHess(coef[measure.par], fMeas, NULL, control=conMeas)
-    oout$hessian[measure.par,measure.par] <- hess1
-  }
-
-  #    vcov <- if (length(coef))
-  #	  solve(oout$hessian)
-  #   else matrix(numeric(0L), 0L, 0L)
-
-  vcov <- matrix(NA, length(coef), length(coef))
-  if (length(coef)) {
-    rrr <- try(solve(oout$hessian), TRUE)
-    if(class(rrr)[1] != "try-error")
-      vcov <- rrr
-  }
-
-  mycoef <- as.list(coef)
-
-  if(!is.CARMA(yuima)){
-    names(mycoef) <- nm
-  }
-  idx.fixed <- orig.idx.fixed
-
-
-
-  mycoef.cont <- mycoef
-  if(length(c(idx.fixed,idx.measure)>0))  # SMI 2/9/14
-    mycoef.cont <- mycoef[-c(idx.fixed,idx.measure)]  # SMI 2/9/14
-
-
-  min.diff <- 0
-  min.jump <- 0
-
-
-  if(length(c(diff.par,drift.par))>0 & !is.CARMA(yuima)){ # LM 04/09/14
-    min.diff <- minusquasilogl(yuima=yuima, param=mycoef[c(diff.par,drift.par)], print=print, env,rcpp=rcpp)
-  }else{
-    if(length(c(diff.par,drift.par))>0 & is.CARMA(yuima)){
[TRUNCATED]

To get the complete diff run:
    svnlook diff /svnroot/yuima -r 745


More information about the Yuima-commits mailing list