[Yuima-commits] r689 - pkg/yuima/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Mar 4 20:20:51 CET 2019
Author: lorenzo
Date: 2019-03-04 20:20:51 +0100 (Mon, 04 Mar 2019)
New Revision: 689
Modified:
pkg/yuima/R/AuxMethodforPPR.R
pkg/yuima/R/qmleLevy.R
Log:
Updated PPR estimation
Modified: pkg/yuima/R/AuxMethodforPPR.R
===================================================================
--- pkg/yuima/R/AuxMethodforPPR.R 2019-03-04 03:54:09 UTC (rev 688)
+++ pkg/yuima/R/AuxMethodforPPR.R 2019-03-04 19:20:51 UTC (rev 689)
@@ -133,8 +133,48 @@
my.envd3 <- new.env()
namesparam<-names(param)
+ resCov<-NULL
if(!(all(namesparam %in% yuimaPPR at PPR@allparamPPR) && length(namesparam)==length(yuimaPPR at PPR@allparamPPR))){
- return(NULL)
+
+ if(length(yuimaPPR at PPR@common)==0){
+ namesCov <-yuimaPPR at PPR@covariates
+ posCov <- length(namesCov)
+ dummydrift <- as.character(yuimaPPR at model@drift[1:posCov])
+ dummydiff0<- NULL
+ for(j in c(1:posCov)){
+ dummydiff0<-c(dummydiff0,
+ as.character(unlist(yuimaPPR at model@diffusion[[j]])))
+ }
+
+ dummydiff <- matrix(dummydiff0, nrow = posCov,
+ ncol = length(dummydiff0)/posCov)
+ dimJump <- length(yuimaPPR at model@jump.coeff[[1]])-length(yuimaPPR at PPR@counting.var)
+ if(dimJump>0){
+ dummyJump0<- NULL
+ for(j in c(1:posCov)){
+ dummyJump0 <- c(dummyJump0,
+ as.character(unlist(yuimaPPR at model@jump.coeff[[j]][1:dimJump])))
+ }
+ dummyJump <- matrix(dummyJump0, nrow=posCov,ncol=dimJump)
+ dummyModel <- setModel(drift = dummydrift,
+ diffusion = dummydiff, jump.coeff =dummyJump,
+ measure = list(df=yuimaPPR at model@measure$df),
+ measure.type = yuimaPPR at model@measure.type[posCov],
+ solve.variable = yuimaPPR at PPR@covariates,
+ state.variable = yuimaPPR at PPR@covariates,
+ xinit=yuimaPPR at model@xinit[posCov])
+ dummydata<-setData(original.data = yuimaPPR at data@original.data[,1:posCov],delta = yuimaPPR at sampling@delta)
+ dummyMod1 <-setYuima(model = dummyModel,
+ data=dummydata)
+ dummyMod1 at sampling<-yuimaPPR at sampling
+
+ resCov <- qmleLevy(yuima = dummyMod1,
+ start=param[dummyMod1 at model@parameter at all],
+ lower = lower[dummyMod1 at model@parameter at all],upper=upper[dummyMod1 at model@parameter at all])
+ }else{
+ return(NULL)
+ }
+ }
}
# construction my.envd1
@@ -379,21 +419,29 @@
fn=Internal.LogLikPPR,
my.envd1=my.envd1,my.envd2=my.envd2,my.envd3=my.envd3),
error=function(){NULL})
+ if(!is.null(Hessian)){
+ Hessian <- Hessian[yuimaPPR at PPR@allparamPPR,yuimaPPR at PPR@allparamPPR]
+ }
cond1 <- my.envd3$YUIMA.PPR at model@solve.variable %in% my.envd3$YUIMA.PPR at PPR@counting.var
cond2 <- diff(as.numeric(my.envd3$YUIMA.PPR at data@original.data[,cond1]))
N.jump <- sum(cond2,na.rm=TRUE)
if(is.null(Hessian)){
- vcov <- matrix(NA,length(out$par),length(out$par))
+ vcov <- matrix(NA,length(out$par[yuimaPPR at PPR@allparamPPR]),
+ length(out$par[yuimaPPR at PPR@allparamPPR]))
}else{
vcov <- solve(Hessian)/N.jump
}
minuslog <- out$value*N.jump
- final_res<-new("yuima.PPR.qmle", call = call, coef = out$par, fullcoef = out$par,
+ final_res<-new("yuima.PPR.qmle", call = call, coef = out$par[yuimaPPR at PPR@allparamPPR],
+ fullcoef = out$par[yuimaPPR at PPR@allparamPPR],
vcov = vcov, min = minuslog, details = out, minuslogl = Internal.LogLikPPR,
method = method, nobs=as.integer(N.jump), model=my.envd3$YUIMA.PPR)
- return(final_res)
-
+ if(!is.null(resCov)){
+ return(list(PPR=final_res,Covariates=resCov))
+ }else{
+ return(final_res)
+ }
}
Modified: pkg/yuima/R/qmleLevy.R
===================================================================
--- pkg/yuima/R/qmleLevy.R 2019-03-04 03:54:09 UTC (rev 688)
+++ pkg/yuima/R/qmleLevy.R 2019-03-04 19:20:51 UTC (rev 689)
@@ -1,258 +1,260 @@
-########################################################################
-# Stepwise estimation for ergodic Levy driven SDE
-########################################################################
-
-
-qmleLevy<-function(yuima,start,lower,upper,joint = FALSE,third = FALSE)
-{
- if(missing(yuima))
- yuima.stop("yuima object is missing.")
-
- if(!is(yuima,"yuima"))
- yuima.stop("This function is for yuima-class.")
-
- sdeModel<-yuima at model
- code <- suppressWarnings(sub("^(.+?)\\(.+", "\\1", sdeModel at measure$df$expr, perl=TRUE))
-
- candinoise<-c("rNIG","rvgamma","rnts","rbgamma")
-
- if(is.na(match(code,candinoise))){
- yuima.stop("This function works only for the standardized normal inverse Gaussian process, variance gamma process, bilateral gamma process, and normal tempered stable process now.")
- }
-
-
- if(length(sdeModel at xinit) == 1){
- args <- unlist(strsplit(suppressWarnings(sub("^.+?\\((.+)\\)", "\\1", sdeModel at measure$df$expr, perl=TRUE)), ","))
- if(code == "rNIG"){
- if(!((abs(eval(parse(text = paste("(",args[5] ,")+(", args[3], ")*(", args[4],
- ")/sqrt((", args[2], ")^2-(", args[3], ")^2)" )))) < 10^(-10))
- && (abs(eval(parse(text = paste("(",args[2], ")^2*(", args[4], ")/(sqrt((", args[2],
- ")^2-(", args[3], ")^2))^3"))) -1) < 10^(-10))))
- {
- yuima.stop("This function is only for standardized Levy noises.")
- }
- }
- else if(code == "rvgamma"){
- if(!((abs(eval(parse(text = paste("(", args[5], ")+2*(", args[2], ")*(", args[4], ")/((", args[3], ")^2-("
- , args[4], ")^2)" )))) < 10^(-10))
- && (abs(eval(parse(text = paste("2*((", args[2], ")*(", args[3], ")^2+(", args[4], ")^2)", "/(("
- , args[3], ")^2-(", args[4], ")^2)^2"))) - 1) < 10^(-10))))
- {
- yuima.stop("This function is only for standardized Levy noises")
- }
- }
- else if((code == "rnts")){
- if(!((abs(eval(parse(text = paste("(", args[6], ")-(", args[2], ")*(",
- args[3], ")*(", args[4], ")^((", args[2],
- ")-1)*gamma(1-(", args[2], "))*(-1/(", args[2],"))*(", args[5],
- ")"
- )))) < 10^(-10))
- &&(abs(eval(parse(text = paste("(", args[3], ")*(", args[2],")*((", args[2], ")-1)*gamma(1-(", args[2], "))*(-1/(", args[2],"))*(", args[4], ")^((",args[2], ")-2)*(", args[5], ")^2-(",
- args[2], ")*(", args[3], ")*(", args[4], ")^((", args[2], ")-1)*gamma(1-(", args[2], "))*(-1/(", args[2],"))"
- ))) - 1) < 10^(-10))))
- {
- yuima.stop("This function is only for standardized Levy processes.")
- }
- }else if(code == "rbgamma"){
- if(!((abs(eval(parse(text = paste("(", args[2], ")/(", args[3],")-(", args[4],")/(", args[5], ")")))) < 10^(-10))
- && (abs(eval(parse(text = paste("(", args[2], ")/(", args[3], ")^2","+(", args[4],")/(", args[5],")^2"
- ))) - 1) < 10^(-10) )))
- {
- yuima.stop("This function is only for standardized Levy processes.")
- }
- }
- }else{
- warning("In this version, the standardized conditions on multidimensional noises can not be verified.
- The expressions of mean and variance are given in help page.")
- # The noise condition checker below does not work now (YU: 3/23).
-
- # args <- suppressWarnings(sub("^.+?\\((.+)\\)", "\\1", sdeModel at measure$df$expr, perl=TRUE))
- # yuimaEnv <- new.env()
- # yuimaEnv$mean <- switch(code,
- # rNIG = function(x=1,alpha,beta,delta0,mu,Lambda){mu+as.vector(delta0/(sqrt(alpha^2-t(beta)%*%Lambda%*%beta)))*Lambda%*%beta},
- # rnts = function(x=1,alpha,a,b,beta,mu,Lambda){mu+gamma(1-alpha)*a*b^(alpha-1)*Lambda%*%beta},
- # rvgamma = function(x=1,lambda,alpha,beta,mu,Lambda){mu+as.vector(2*lambda/(alpha^2-t(beta)%*%Lambda%*%beta)^2)*beta}
- # )
- #
- # yuimaEnv$covariance <- switch(code,
- # rNIG = function(x=1,alpha,beta,delta0,mu,Lambda){as.vector(delta0/(sqrt(alpha^2-t(beta)%*%Lambda%*%beta))^3)*Lambda%*%beta%*%t(beta)%*%Lambda+as.vector(delta0/sqrt(alpha^2-t(beta)%*%Lambda%*%beta))*Lambda},
- # rnts = function(x=1,alpha,a,b,beta,mu,Lambda){a*(1-alpha)*gamma(1-alpha)*b^(alpha-2)*Lambda%*%beta%*%t(beta)%*%Lambda+a*gamma(1-alpha)*b^(alpha-1)*Lambda},
- # rvgamma = function(x=1,lambda,alpha,beta,mu,Lambda){as.vector(4*lambda/(alpha^2-t(beta)%*%Lambda%*%beta)^2)*Lambda%*%beta%*%t(beta)%*%Lambda+as.vector(2*lambda/alpha^2-t(beta)%*%Lambda%*%beta)*Lambda}
- # )
- # judgemean<-sum(eval(parse(text = paste("mean","(",args,")")),yuimaEnv)==numeric(length(sdeModel at xinit)))
- # judgecovariance<-sum(eval(eval(parse(text = paste("covariance","(",args,")")),yuimaEnv)==diag(1,length(sdeModel at xinit))))
- # if(!((judgemean==length(sdeModel at xinit))&&(judgecovariance==length(sdeModel at xinit)*length(sdeModel at xinit))))
- # {
- # yuima.stop("This function is only for standardized Levy processes.")
- # }
- }
- yuima at sampling@delta <- yuima at sampling@delta[1]
- yuima at model@noise.number <- as.integer(yuima at model@equation.number)
- if(!joint){
- DRIFT <- yuima at model@drift
- DRPAR <- yuima at model@parameter at drift
-
- if(length(yuima at model@parameter at jump)>0)
- fullcoef <- yuima at model@parameter at jump
-
- if(length(DRPAR)>0)
- fullcoef <- c(fullcoef, DRPAR)
-
- oo <- match(yuima at model@parameter at all, fullcoef)
- yuima at model@parameter at all <- yuima at model@parameter at all[order(oo)]
-
- oo <- match(names(start), fullcoef)
- start <- start[order(oo)]
-
- oo <- match(names(upper), fullcoef)
- upper <- upper[order(oo)]
-
- oo <- match(names(lower), fullcoef)
- lower <- lower[order(oo)]
-
- yuima at model@diffusion <- yuima at model@jump.coeff
- yuima at model@parameter at diffusion <- yuima at model@parameter at jump[1:length(yuima at model@parameter at jump)]
- yuima at model@parameter at all <- yuima at model@parameter at diffusion
-
- for(i in 1:length(yuima at model@drift)){
- yuima at model@drift[i] <- expression((0))
- }
-
- yuima at model@jump.coeff <- list()
- yuima at model@parameter at drift <- character(0)
- yuima at model@measure <- list()
- yuima at model@jump.variable <- character(0)
- yuima at model@measure.type <- character(0)
- yuima at model@parameter at jump <- character(0)
- yuima at model@parameter at measure <- character(0)
-
- diffstart <- start[1:length(yuima at model@parameter at diffusion)]
- diffupper <- upper[1:length(yuima at model@parameter at diffusion)]
- difflower <- lower[1:length(yuima at model@parameter at diffusion)]
- fres <- qmle(yuima=yuima,start=diffstart,lower=difflower,upper=diffupper,rcpp = TRUE,joint = FALSE,method = "L-BFGS-B")
-
- DIPAR <- yuima at model@parameter at diffusion
- DIFFUSION <- yuima at model@diffusion
- yuima at model@parameter at all <- DRPAR
- yuima at model@parameter at drift <- DRPAR
- yuima at model@drift <- DRIFT
-
- dristart <- start[-(1:length(yuima at model@parameter at diffusion))]
- driupper <- upper[-(1:length(yuima at model@parameter at diffusion))]
- drilower <- lower[-(1:length(yuima at model@parameter at diffusion))]
-
- partcoef <- yuima at model@parameter at diffusion
- ov <- match(yuima at model@parameter at drift,partcoef)
- ovp <- which(!is.na(ov))
- if(length(ovp)>0)
- {yuima at model@parameter at drift <- yuima at model@parameter at drift[-ovp]}
- yuima at model@parameter at all <- yuima at model@parameter at drift
-
- ma <- match(names(fres at coef),partcoef)
- sort <- fres at coef[order(ma)]
- esti <- numeric(length(partcoef))
- newdiff <- yuima at model@diffusion
- newdri <- yuima at model@drift
-
- for(i in 1:length(partcoef))
- {
- esti[i] <- as.character(fres at coef[[i]])
- }
-
- if(length(yuima at model@drift) == 1){
- for(i in 1:length(partcoef))
- {
- newdri <- gsub(partcoef[i],esti[i],newdri)
- yuima at model@drift[1] <- parse(text = newdri[1])
- newdiff[[1]] <- gsub(partcoef[i],esti[i],newdiff[[1]])
- yuima at model@diffusion[[1]] <- parse(text = newdiff[[1]])
- }
- }else{
- for(i in 1:length(partcoef))
- {
- for(j in 1:length(yuima at model@drift))
- {
- newdri[j] <- gsub(partcoef[i],esti[i],newdri[j])
- yuima at model@drift[j] <- parse(text = newdri[j])
- }
- for(k in 1:length(yuima at model@diffusion)){
- for(l in 1:length(yuima at model@diffusion[[1]])){
- newdiff[[k]][l] <- gsub(partcoef[i],esti[i],newdiff[[k]][l])
- yuima at model@diffusion[[k]][l] <- parse(text = newdiff[[k]][l])
- }
- }
- }
- }
- yuima at model@parameter at diffusion <- character(0)
-
- sres<-qmle(yuima=yuima,start=dristart,lower=drilower,upper=driupper,rcpp = TRUE,method = "L-BFGS-B")
-
- if((length(ovp) == 0) && (third)){
- yuima at model@diffusion <- DIFFUSION
- yuima at model@drift <- DRIFT
- yuima at model@parameter at diffusion <- DIPAR
- yuima at model@parameter at all <- DIPAR
- newdri <- yuima at model@drift
-
- for(i in 1:length(sres at coef))
- {
- esti[i] <- as.character(sres at coef[[i]])
- }
-
- if(length(yuima at model@drift)==1){
- for(i in 1:length(sres at coef))
- {
- newdri <- gsub(yuima at model@parameter at drift[i],esti[i],newdri)
- yuima at model@drift[1] <- parse(text=newdri[1])
- }
- }else{
- for(i in 1:length(sres at coef))
- {
- for(j in 1:length(yuima at model@drift))
- {
- newdri[j] <- gsub(yuima at model@parameter at drift[i],esti[i],newdri[j])
- yuima at model@drift[j] <- parse(text = newdri[j])
- }
- }
- }
-
- yuima at model@parameter at drift <- character(0)
- too <- match(names(fres at coef),names(diffstart))
- diffstart <- diffstart[order(too)]
- for(i in 1:length(diffstart)){
- diffstart[[i]] <- fres at coef[[i]]
- }
- tres <- qmle(yuima=yuima,start=diffstart,lower=difflower,upper=diffupper,rcpp = TRUE,joint = FALSE,method = "L-BFGS-B")
-
- res <- list(first = fres at coef, second = sres at coef, third = tres at coef)
- }else if((length(ovp) > 0) || !(third)){
- res <- list(first = fres at coef, second = sres at coef)}
- else{
- yuima.stop("third estimation may be theoretical invalid under the presence of an overlapping parameter.")
- }
- }else{
- if(third){
- yuima.stop("third estimation does not make sense in joint estimation.")
- }
- if(length(yuima at model@parameter at jump)>0)
- fullcoef <- yuima at model@parameter at jump
- if(length(yuima at model@parameter at drift)>0)
- fullcoef <- c(fullcoef, yuima at model@parameter at drift)
- oo <- match(yuima at model@parameter at all, fullcoef)
- yuima at model@parameter at all <- yuima at model@parameter at all[order(oo)]
- yuima at model@parameter at all <- yuima at model@parameter at all[1:length(which(!is.na(oo)))]
- yuima at model@diffusion <- yuima at model@jump.coeff
- yuima at model@jump.coeff <- list()
- yuima at model@parameter at diffusion <- yuima at model@parameter at jump[1:length(yuima at model@parameter at jump)]
- yuima at model@measure <- list()
- yuima at model@jump.variable <- character(0)
- yuima at model@measure.type <- character(0)
- yuima at model@parameter at jump <- character(0)
- yuima at model@parameter at measure <- character(0)
- jres<-qmle(yuima,start = start,lower = lower,upper = upper,rcpp = TRUE, joint = TRUE,method = "L-BFGS-B")
- res<-list(joint = jres at coef)
- }
- res
-}
-
-
+########################################################################
+# Stepwise estimation for ergodic Levy driven SDE
+########################################################################
+
+
+qmleLevy<-function(yuima,start,lower,upper,joint = FALSE,third = FALSE)
+{
+ if(missing(yuima))
+ yuima.stop("yuima object is missing.")
+
+ if(!is(yuima,"yuima"))
+ yuima.stop("This function is for yuima-class.")
+
+ sdeModel<-yuima at model
+ if(class(sdeModel at measure$df)!="yuima.law"){
+ code <- suppressWarnings(sub("^(.+?)\\(.+", "\\1", sdeModel at measure$df$expr, perl=TRUE))
+
+ candinoise<-c("rNIG","rvgamma","rnts","rbgamma")
+
+ if(is.na(match(code,candinoise))){
+ yuima.stop("This function works only for the standardized normal inverse Gaussian process, variance gamma process, bilateral gamma process, and normal tempered stable process now.")
+ }
+
+
+ if(length(sdeModel at xinit) == 1){
+ args <- unlist(strsplit(suppressWarnings(sub("^.+?\\((.+)\\)", "\\1", sdeModel at measure$df$expr, perl=TRUE)), ","))
+ if(code == "rNIG"){
+ if(!((abs(eval(parse(text = paste("(",args[5] ,")+(", args[3], ")*(", args[4],
+ ")/sqrt((", args[2], ")^2-(", args[3], ")^2)" )))) < 10^(-10))
+ && (abs(eval(parse(text = paste("(",args[2], ")^2*(", args[4], ")/(sqrt((", args[2],
+ ")^2-(", args[3], ")^2))^3"))) -1) < 10^(-10))))
+ {
+ yuima.stop("This function is only for standardized Levy noises.")
+ }
+ }
+ else if(code == "rvgamma"){
+ if(!((abs(eval(parse(text = paste("(", args[5], ")+2*(", args[2], ")*(", args[4], ")/((", args[3], ")^2-("
+ , args[4], ")^2)" )))) < 10^(-10))
+ && (abs(eval(parse(text = paste("2*((", args[2], ")*(", args[3], ")^2+(", args[4], ")^2)", "/(("
+ , args[3], ")^2-(", args[4], ")^2)^2"))) - 1) < 10^(-10))))
+ {
+ yuima.stop("This function is only for standardized Levy noises")
+ }
+ }
+ else if((code == "rnts")){
+ if(!((abs(eval(parse(text = paste("(", args[6], ")-(", args[2], ")*(",
+ args[3], ")*(", args[4], ")^((", args[2],
+ ")-1)*gamma(1-(", args[2], "))*(-1/(", args[2],"))*(", args[5],
+ ")"
+ )))) < 10^(-10))
+ &&(abs(eval(parse(text = paste("(", args[3], ")*(", args[2],")*((", args[2], ")-1)*gamma(1-(", args[2], "))*(-1/(", args[2],"))*(", args[4], ")^((",args[2], ")-2)*(", args[5], ")^2-(",
+ args[2], ")*(", args[3], ")*(", args[4], ")^((", args[2], ")-1)*gamma(1-(", args[2], "))*(-1/(", args[2],"))"
+ ))) - 1) < 10^(-10))))
+ {
+ yuima.stop("This function is only for standardized Levy processes.")
+ }
+ }else if(code == "rbgamma"){
+ if(!((abs(eval(parse(text = paste("(", args[2], ")/(", args[3],")-(", args[4],")/(", args[5], ")")))) < 10^(-10))
+ && (abs(eval(parse(text = paste("(", args[2], ")/(", args[3], ")^2","+(", args[4],")/(", args[5],")^2"
+ ))) - 1) < 10^(-10) )))
+ {
+ yuima.stop("This function is only for standardized Levy processes.")
+ }
+ }
+ }else{
+ warning("In this version, the standardized conditions on multidimensional noises can not be verified.
+ The expressions of mean and variance are given in help page.")
+ # The noise condition checker below does not work now (YU: 3/23).
+
+ # args <- suppressWarnings(sub("^.+?\\((.+)\\)", "\\1", sdeModel at measure$df$expr, perl=TRUE))
+ # yuimaEnv <- new.env()
+ # yuimaEnv$mean <- switch(code,
+ # rNIG = function(x=1,alpha,beta,delta0,mu,Lambda){mu+as.vector(delta0/(sqrt(alpha^2-t(beta)%*%Lambda%*%beta)))*Lambda%*%beta},
+ # rnts = function(x=1,alpha,a,b,beta,mu,Lambda){mu+gamma(1-alpha)*a*b^(alpha-1)*Lambda%*%beta},
+ # rvgamma = function(x=1,lambda,alpha,beta,mu,Lambda){mu+as.vector(2*lambda/(alpha^2-t(beta)%*%Lambda%*%beta)^2)*beta}
+ # )
+ #
+ # yuimaEnv$covariance <- switch(code,
+ # rNIG = function(x=1,alpha,beta,delta0,mu,Lambda){as.vector(delta0/(sqrt(alpha^2-t(beta)%*%Lambda%*%beta))^3)*Lambda%*%beta%*%t(beta)%*%Lambda+as.vector(delta0/sqrt(alpha^2-t(beta)%*%Lambda%*%beta))*Lambda},
+ # rnts = function(x=1,alpha,a,b,beta,mu,Lambda){a*(1-alpha)*gamma(1-alpha)*b^(alpha-2)*Lambda%*%beta%*%t(beta)%*%Lambda+a*gamma(1-alpha)*b^(alpha-1)*Lambda},
+ # rvgamma = function(x=1,lambda,alpha,beta,mu,Lambda){as.vector(4*lambda/(alpha^2-t(beta)%*%Lambda%*%beta)^2)*Lambda%*%beta%*%t(beta)%*%Lambda+as.vector(2*lambda/alpha^2-t(beta)%*%Lambda%*%beta)*Lambda}
+ # )
+ # judgemean<-sum(eval(parse(text = paste("mean","(",args,")")),yuimaEnv)==numeric(length(sdeModel at xinit)))
+ # judgecovariance<-sum(eval(eval(parse(text = paste("covariance","(",args,")")),yuimaEnv)==diag(1,length(sdeModel at xinit))))
+ # if(!((judgemean==length(sdeModel at xinit))&&(judgecovariance==length(sdeModel at xinit)*length(sdeModel at xinit))))
+ # {
+ # yuima.stop("This function is only for standardized Levy processes.")
+ # }
+ }
+ }else{fullcoef<-NULL}
+ yuima at sampling@delta <- yuima at sampling@delta[1]
+ yuima at model@noise.number <- as.integer(yuima at model@equation.number)
+ if(!joint){
+ DRIFT <- yuima at model@drift
+ DRPAR <- yuima at model@parameter at drift
+
+ if(length(yuima at model@parameter at jump)>0)
+ fullcoef <- yuima at model@parameter at jump
+
+ if(length(DRPAR)>0)
+ fullcoef <- c(fullcoef, DRPAR)
+
+ oo <- match(yuima at model@parameter at all, fullcoef)
+ yuima at model@parameter at all <- yuima at model@parameter at all[order(oo)]
+
+ oo <- match(names(start), fullcoef)
+ start <- start[order(oo)]
+
+ oo <- match(names(upper), fullcoef)
+ upper <- upper[order(oo)]
+
+ oo <- match(names(lower), fullcoef)
+ lower <- lower[order(oo)]
+
+ yuima at model@diffusion <- yuima at model@jump.coeff
+ yuima at model@parameter at diffusion <- yuima at model@parameter at jump[1:length(yuima at model@parameter at jump)]
+ yuima at model@parameter at all <- yuima at model@parameter at diffusion
+
+ for(i in 1:length(yuima at model@drift)){
+ yuima at model@drift[i] <- expression((0))
+ }
+
+ yuima at model@jump.coeff <- list()
+ yuima at model@parameter at drift <- character(0)
+ yuima at model@measure <- list()
+ yuima at model@jump.variable <- character(0)
+ yuima at model@measure.type <- character(0)
+ yuima at model@parameter at jump <- character(0)
+ yuima at model@parameter at measure <- character(0)
+
+ diffstart <- start[1:length(yuima at model@parameter at diffusion)]
+ diffupper <- upper[1:length(yuima at model@parameter at diffusion)]
+ difflower <- lower[1:length(yuima at model@parameter at diffusion)]
+ fres <- qmle(yuima=yuima,start=diffstart,lower=difflower,upper=diffupper,rcpp = TRUE,joint = FALSE,method = "L-BFGS-B")
+
+ DIPAR <- yuima at model@parameter at diffusion
+ DIFFUSION <- yuima at model@diffusion
+ yuima at model@parameter at all <- DRPAR
+ yuima at model@parameter at drift <- DRPAR
+ yuima at model@drift <- DRIFT
+
+ dristart <- start[-(1:length(yuima at model@parameter at diffusion))]
+ driupper <- upper[-(1:length(yuima at model@parameter at diffusion))]
+ drilower <- lower[-(1:length(yuima at model@parameter at diffusion))]
+
+ partcoef <- yuima at model@parameter at diffusion
+ ov <- match(yuima at model@parameter at drift,partcoef)
+ ovp <- which(!is.na(ov))
+ if(length(ovp)>0)
+ {yuima at model@parameter at drift <- yuima at model@parameter at drift[-ovp]}
+ yuima at model@parameter at all <- yuima at model@parameter at drift
+
+ ma <- match(names(fres at coef),partcoef)
+ sort <- fres at coef[order(ma)]
+ esti <- numeric(length(partcoef))
+ newdiff <- yuima at model@diffusion
+ newdri <- yuima at model@drift
+
+ for(i in 1:length(partcoef))
+ {
+ esti[i] <- as.character(fres at coef[[i]])
+ }
+
+ if(length(yuima at model@drift) == 1){
+ for(i in 1:length(partcoef))
+ {
+ newdri <- gsub(partcoef[i],esti[i],newdri)
+ yuima at model@drift[1] <- parse(text = newdri[1])
+ newdiff[[1]] <- gsub(partcoef[i],esti[i],newdiff[[1]])
+ yuima at model@diffusion[[1]] <- parse(text = newdiff[[1]])
+ }
+ }else{
+ for(i in 1:length(partcoef))
+ {
+ for(j in 1:length(yuima at model@drift))
+ {
+ newdri[j] <- gsub(partcoef[i],esti[i],newdri[j])
+ yuima at model@drift[j] <- parse(text = newdri[j])
+ }
+ for(k in 1:length(yuima at model@diffusion)){
+ for(l in 1:length(yuima at model@diffusion[[1]])){
+ newdiff[[k]][l] <- gsub(partcoef[i],esti[i],newdiff[[k]][l])
+ yuima at model@diffusion[[k]][l] <- parse(text = newdiff[[k]][l])
+ }
+ }
+ }
+ }
+ yuima at model@parameter at diffusion <- character(0)
+
+ sres<-qmle(yuima=yuima,start=dristart,lower=drilower,upper=driupper,rcpp = TRUE,method = "L-BFGS-B")
+
+ if((length(ovp) == 0) && (third)){
+ yuima at model@diffusion <- DIFFUSION
+ yuima at model@drift <- DRIFT
+ yuima at model@parameter at diffusion <- DIPAR
+ yuima at model@parameter at all <- DIPAR
+ newdri <- yuima at model@drift
+
+ for(i in 1:length(sres at coef))
+ {
+ esti[i] <- as.character(sres at coef[[i]])
+ }
+
+ if(length(yuima at model@drift)==1){
+ for(i in 1:length(sres at coef))
+ {
+ newdri <- gsub(yuima at model@parameter at drift[i],esti[i],newdri)
+ yuima at model@drift[1] <- parse(text=newdri[1])
+ }
+ }else{
+ for(i in 1:length(sres at coef))
+ {
+ for(j in 1:length(yuima at model@drift))
+ {
+ newdri[j] <- gsub(yuima at model@parameter at drift[i],esti[i],newdri[j])
+ yuima at model@drift[j] <- parse(text = newdri[j])
+ }
+ }
+ }
+
+ yuima at model@parameter at drift <- character(0)
+ too <- match(names(fres at coef),names(diffstart))
+ diffstart <- diffstart[order(too)]
+ for(i in 1:length(diffstart)){
+ diffstart[[i]] <- fres at coef[[i]]
+ }
+ tres <- qmle(yuima=yuima,start=diffstart,lower=difflower,upper=diffupper,rcpp = TRUE,joint = FALSE,method = "L-BFGS-B")
+
+ res <- list(first = fres at coef, second = sres at coef, third = tres at coef)
+ }else if((length(ovp) > 0) || !(third)){
+ res <- list(first = fres at coef, second = sres at coef)}
+ else{
+ yuima.stop("third estimation may be theoretical invalid under the presence of an overlapping parameter.")
+ }
+ }else{
+ if(third){
+ yuima.stop("third estimation does not make sense in joint estimation.")
+ }
+ if(length(yuima at model@parameter at jump)>0)
+ fullcoef <- yuima at model@parameter at jump
+ if(length(yuima at model@parameter at drift)>0)
+ fullcoef <- c(fullcoef, yuima at model@parameter at drift)
+ oo <- match(yuima at model@parameter at all, fullcoef)
+ yuima at model@parameter at all <- yuima at model@parameter at all[order(oo)]
+ yuima at model@parameter at all <- yuima at model@parameter at all[1:length(which(!is.na(oo)))]
+ yuima at model@diffusion <- yuima at model@jump.coeff
+ yuima at model@jump.coeff <- list()
+ yuima at model@parameter at diffusion <- yuima at model@parameter at jump[1:length(yuima at model@parameter at jump)]
+ yuima at model@measure <- list()
+ yuima at model@jump.variable <- character(0)
+ yuima at model@measure.type <- character(0)
+ yuima at model@parameter at jump <- character(0)
+ yuima at model@parameter at measure <- character(0)
+ jres<-qmle(yuima,start = start,lower = lower,upper = upper,rcpp = TRUE, joint = TRUE,method = "L-BFGS-B")
+ res<-list(joint = jres at coef)
+ }
+ res
+}
+
+
More information about the Yuima-commits
mailing list