[Yuima-commits] r596 - pkg/yuima/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Mar 27 07:54:14 CEST 2017
Author: eguchi
Date: 2017-03-27 07:54:13 +0200 (Mon, 27 Mar 2017)
New Revision: 596
Added:
pkg/yuima/R/qmleLevy.R
Log:
added
Added: pkg/yuima/R/qmleLevy.R
===================================================================
--- pkg/yuima/R/qmleLevy.R (rev 0)
+++ pkg/yuima/R/qmleLevy.R 2017-03-27 05:54:13 UTC (rev 596)
@@ -0,0 +1,258 @@
+########################################################################
+# 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(!((round(eval(parse(text = paste("(",args[5] ,")+(", args[3], ")*(", args[4],
+ ")/sqrt((", args[2], ")^2-(", args[3], ")^2)" ))), digit = 10) == 0)
+ && (round(eval(parse(text = paste("(",args[2], ")^2*(", args[4], ")/(sqrt((", args[2],
+ ")^2-(", args[3], ")^2))^3"))), digit = 10) == 1)))
+ {
+ yuima.stop("This function is only for standardized Levy noises.")
+ }
+ }
+ else if(code == "rvgamma"){
+ if(!((round(eval(parse(text = paste("(", args[5], ")+2*(", args[2], ")*(", args[4], ")/((", args[3], ")^2-("
+ , args[4], ")^2)" ))), digit = 10) == 0)
+ && (round(eval(parse(text = paste("2*((", args[2], ")*(", args[3], ")^2+(", args[4], ")^2)", "/(("
+ , args[3], ")^2-(", args[4], ")^2)^2"))), digit = 10) == 1)))
+ {
+ yuima.stop("This function is only for standardized Levy noises")
+ }
+ }
+ else if((code == "rnts")){
+ if(!((round(eval(parse(text = paste("(", args[6], ")-(", args[2], ")*(",
+ args[3], ")*(", args[4], ")^((", args[2],
+ ")-1)*gamma(1-(", args[2], "))*(-1/(", args[2],"))*(", args[5],
+ ")"
+ ))), digit = 10) == 0)
+ &&(round(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],"))"
+ ))), digit = 10) == 1)))
+ {
+ yuima.stop("This function is only for standardized Levy processes.")
+ }
+ }else if(code == "rbgamma"){
+ if(!((round(eval(parse(text = paste("(", args[2], ")/(", args[3],")-(", args[4],")/(", args[5], ")"))),digit = 10) == 0)
+ && (round(eval(parse(text = paste("(", args[2], ")/(", args[3], ")^2","+(", args[4],")/(", args[5],")^2"
+ ))), digit = 10) == 1)))
+ {
+ 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
+}
+
+
More information about the Yuima-commits
mailing list