[Yuima-commits] r501 - pkg/yuima/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Oct 29 23:07:13 CEST 2016
Author: lorenzo
Date: 2016-10-29 23:07:13 +0200 (Sat, 29 Oct 2016)
New Revision: 501
Modified:
pkg/yuima/R/qmle.R
pkg/yuima/R/simulateForPpr.R
Log:
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2016-10-29 21:03:26 UTC (rev 500)
+++ pkg/yuima/R/qmle.R 2016-10-29 21:07:13 UTC (rev 501)
@@ -123,6 +123,15 @@
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")
@@ -1402,7 +1411,7 @@
# cov[unique(c(measure.par)),unique(c(measure.par))]<-dummycovCarmaNoise
# carma_final_res<-list(mle=final_res,Incr=inc.levy,model=yuima)
- if(Est.Incr=="Carma.IncPar"){
+ if(Est.Incr=="Carma.IncPar"||Est.Incr=="IncrPar"){
#inc.levy.fin<-zoo(inc.levy,tt,frequency=1/env$h)
inc.levy.fin<-zoo(inc.levy,tt[(1+length(tt)-length(inc.levy)):length(tt)])
carma_final_res<-new("yuima.carma.qmle", call = call, coef = coef, fullcoef = unlist(coef),
@@ -1411,7 +1420,7 @@
model = yuima at model, nobs=yuima.nobs,
logL.Incr = tryCatch(-result.Lev$value,error=function(theta){NULL}))
}else{
- if(Est.Incr=="Carma.Par"){
+ if(Est.Incr=="Carma.Par"||Est.Incr=="NoIncr"){
carma_final_res<-new("mle", call = call, coef = coef, fullcoef = unlist(coef),
vcov = cov, min = min, details = oout, minuslogl = minusquasilogl,
method = method, nobs=yuima.nobs)
Modified: pkg/yuima/R/simulateForPpr.R
===================================================================
--- pkg/yuima/R/simulateForPpr.R 2016-10-29 21:03:26 UTC (rev 500)
+++ pkg/yuima/R/simulateForPpr.R 2016-10-29 21:07:13 UTC (rev 501)
@@ -51,6 +51,60 @@
increment.L = increment.L, method = method, hurst = hurst,
methodfGn = methodfGn, sampling = sampling,
subsampling = subsampling){
+ ROLDVER<-TRUE
+ if(ROLDVER){
+ object <- aux.simulatPprROldVersion(object, nsim = nsim, seed = seed,
+ xinit = xinit, true.parameter = true.parameter,
+ space.discretized = space.discretized, increment.W = increment.W,
+ increment.L = increment.L, method = method, hurst = hurst,
+ methodfGn = methodfGn, sampling = sampling,
+ subsampling = subsampling)
+ }else{
+ object <- aux.simulatPprROldNew(object, nsim = nsim, seed = seed,
+ xinit = xinit, true.parameter = true.parameter,
+ space.discretized = space.discretized, increment.W = increment.W,
+ increment.L = increment.L, method = method, hurst = hurst,
+ methodfGn = methodfGn, sampling = sampling,
+ subsampling = subsampling)
+ }
+ return(object)
+}
+
+aux.simulatPprROldNew<-function(object, nsim = nsim, seed = seed,
+ xinit = xinit, true.parameter = true.parameter,
+ space.discretized = space.discretized, increment.W = increment.W,
+ increment.L = increment.L, method = method, hurst = hurst,
+ methodfGn = methodfGn, sampling = sampling,
+ subsampling = subsampling){
+
+ # We need an expression for the evaluation of the hazard
+ is(object)
+ yuimaPpr<-object
+
+ envPpr <- list()
+
+ dY <- paste0("d",yuimaPpr at Ppr@var.dx)
+ Time <- sampling at grid[[1]]
+
+ IntKernExpr<- function(Kern,dy){
+ dum<-paste(Kern,dy,sep="*")
+ dum<-paste0(dum, collapse = " + ")
+ dum <- paste0("sum( ", dum, " )")
+ return(parse(text = dum))
+ }
+ IntegKern <- lapply(yuimaPpr at Kernel@Integrand at IntegrandList,IntKernExpr,dY)
+
+
+ return(NULL)
+}
+
+
+aux.simulatPprROldVersion <- function(object, nsim = nsim, seed = seed,
+ xinit = xinit, true.parameter = true.parameter,
+ space.discretized = space.discretized, increment.W = increment.W,
+ increment.L = increment.L, method = method, hurst = hurst,
+ methodfGn = methodfGn, sampling = sampling,
+ subsampling = subsampling){
Time <- sampling at Terminal
numbVardx <- length(object at Ppr@var.dx)
numbCountVar <- length(object at Ppr@counting.var)
@@ -70,60 +124,60 @@
assign("t",object at gFun@param at time.var, envir = my.env)
- nameu <- object at gFun@param at time.var
- assign("dt",sampling at delta, envir = my.env)
+ nameu <- object at gFun@param at time.var
+ assign("dt",sampling at delta, envir = my.env)
- if(is.null(increment.W)){
- dimW <- length(object at model@diffusion[[1]])
- W <- matrix(rnorm(dimW*sampling at n,mean=0,sd= sqrt(sampling at delta)),nrow=dimW,ncol=sampling at n)
- }
- Condcovariate <- TRUE
- if(is.null(increment.L)){
- dimL <- length(object at model@jump.coeff[[1]])
- L <- matrix(0,nrow=dimL,ncol=sampling at n)
- Condcovariate <- FALSE
- # if(length(object at Ppr@covariates)!=0)
- # Condcovariate <- TRUE
- cond <- !(object at model@solve.variable %in% object at Ppr@counting.var)
- if(any(cond)){
- Condcovariate <- TRUE
- }
- dimMd <- length(object at model@solve.variable)
- dumMod <- setModel(drift = rep("0",dimMd),
- diffusion = matrix("0",dimMd,1),
- jump.coeff = diag("1",dimMd,dimMd),
- measure = object at Ppr@Info.measure$measure,
- measure.type = object at Ppr@Info.measure$type,
- solve.variable = object at model@solve.variable)
- if(length(object at model@parameter at measure)!=0){
- simMod <- simulate(object = dumMod,
- true.parameter = true.parameter[object at model@parameter at measure],
- sampling = sampling)
- }else{
- simMod <- simulate(object = dumMod,
- sampling = sampling)
- }
+ if(is.null(increment.W)){
+ dimW <- length(object at model@diffusion[[1]])
+ W <- matrix(rnorm(dimW*sampling at n,mean=0,sd= sqrt(sampling at delta)),nrow=dimW,ncol=sampling at n)
+ }
+ Condcovariate <- TRUE
+ if(is.null(increment.L)){
+ dimL <- length(object at model@jump.coeff[[1]])
+ L <- matrix(0,nrow=dimL,ncol=sampling at n)
+ Condcovariate <- FALSE
+ # if(length(object at Ppr@covariates)!=0)
+ # Condcovariate <- TRUE
+ cond <- !(object at model@solve.variable %in% object at Ppr@counting.var)
+ if(any(cond)){
+ Condcovariate <- TRUE
+ }
+ dimMd <- length(object at model@solve.variable)
+ dumMod <- setModel(drift = rep("0",dimMd),
+ diffusion = matrix("0",dimMd,1),
+ jump.coeff = diag("1",dimMd,dimMd),
+ measure = object at Ppr@Info.measure$measure,
+ measure.type = object at Ppr@Info.measure$type,
+ solve.variable = object at model@solve.variable)
+ if(length(object at model@parameter at measure)!=0){
+ simMod <- simulate(object = dumMod,
+ true.parameter = true.parameter[object at model@parameter at measure],
+ sampling = sampling)
+ }else{
+ simMod <- simulate(object = dumMod,
+ sampling = sampling)
+ }
- L <- t(diff(simMod at data@original.data))
- }
+ L <- t(diff(simMod at data@original.data))
+ }
assign("Condcovariate",Condcovariate, envir = my.env)
assign("W", W, envir = my.env)
- rownames(L)<- object at model@solve.variable
+ rownames(L)<- object at model@solve.variable
- assign("L", L, envir = my.env)
+ assign("L", L, envir = my.env)
- assign("All.labKern",object at Kernel@variable.Integral,envir = my.env)
- assign("All.labgFun",object at gFun@param,envir = my.env)
+ assign("All.labKern",object at Kernel@variable.Integral,envir = my.env)
+ assign("All.labgFun",object at gFun@param,envir = my.env)
- Fun1 <- function(u,env){
- part <- seq(0,u,by=env$dt)
- env$t<-part[-length(part)]
- if(Condcovariate){
+ Fun1 <- function(u,env){
+ part <- seq(0,u,by=env$dt)
+ env$t<-part[-length(part)]
+ if(Condcovariate){
yuima<- object at model
for(i in c(1:length(object at Ppr@covariates))){
assign(object at Ppr@covariates[i],
@@ -140,8 +194,8 @@
# sum(yuima at solve.variable!=object at Ppr@covariates), dim(Linc)[2])
Linc[yuima at solve.variable!=object at Ppr@covariates,] <- 0
DumUnderlMod <- simulate(yuima, true.parameter = true.parameter,
- increment.L = env$L[,c(1:(length(part)-1))],
- sampling = setSampling(Terminal = u, n= (length(part)-1)))
+ increment.L = env$L[,c(1:(length(part)-1))],
+ sampling = setSampling(Terminal = u, n= (length(part)-1)))
for(i in c(1:length(object at Ppr@covariates))){
@@ -149,77 +203,77 @@
assign(object at Ppr@covariates[i], as.numeric(VariableDum), envir = env)
}
}
- }
- (log(env$U)+sum(eval(env$gFun,envir = env)*env$dt))^2
- }
+ }
+ (log(env$U)+sum(eval(env$gFun,envir = env)*env$dt))^2
+ }
- Fun2 <- function(u,env){
- u <- max(env$old_u,u)
- dumpart <- seq(0,env$old_u, by=env$dt)
- part <- seq(env$old_u,u,by=env$dt)
- t_k <- env$t
- env$t<-part[-length(part)]
- if(u>=sampling at Terminal){
- # Think a better solution
- my.env$utrue<-u
- return(0)
- }
- if(Condcovariate){
- LevIncr <- env$L[, length(dumpart)+c(1:(length(env$t)))]
- LevIncr[object at Ppr@counting.var,]<-0
- yuima<- object at model
- xinit<- numeric(length(object at Ppr@covariates))
- names(xinit)<- object at Ppr@covariates
- for(i in c(1:length(object at Ppr@covariates))){
+ Fun2 <- function(u,env){
+ u <- max(env$old_u,u)
+ dumpart <- seq(0,env$old_u, by=env$dt)
+ part <- seq(env$old_u,u,by=env$dt)
+ t_k <- env$t
+ env$t<-part[-length(part)]
+ if(u>=sampling at Terminal){
+ # Think a better solution
+ my.env$utrue<-u
+ return(0)
+ }
+ if(Condcovariate){
+ LevIncr <- env$L[, length(dumpart)+c(1:(length(env$t)))]
+ LevIncr[object at Ppr@counting.var,]<-0
+ yuima<- object at model
+ xinit<- numeric(length(object at Ppr@covariates))
+ names(xinit)<- object at Ppr@covariates
+ for(i in c(1:length(object at Ppr@covariates))){
xinit[i] <- env[[object at Ppr@covariates[i]]]
- }
+ }
- xinitCount <- numeric(length(object at Ppr@counting.var))
- names(xinitCount) <- object at Ppr@counting.var
- for(i in c(1:length(xinitCount))){
- xinitCount[i] <- tail(env[[object at Ppr@counting.var[i]]],n = 1)
- }
- xinit <- c(xinit,xinitCount)
- if(part[length(part)]-part[1]!=0){
- DumVarCov <- simulate(yuima,
- true.parameter = true.parameter,
- increment.L = LevIncr,
- sampling = setSampling(Terminal = (part[length(part)]-part[1]),
- n = dim(LevIncr)[2]),
- xinit=xinit[yuima at solve.variable])
- for(i in c(1:length(object at Ppr@covariates))){
- VariableDum <- DumVarCov at data@original.data[,yuima at solve.variable==object at Ppr@covariates[i]]
- assign(object at Ppr@covariates[i], as.numeric(VariableDum), envir = env)
- }
- }else{
- for(i in c(1:length(object at Ppr@covariates))){
- VariableDum <- xinit[yuima at solve.variable==object at Ppr@covariates[i]]
- assign(object at Ppr@covariates[i], as.numeric(VariableDum), envir = env)
- }
- }
- #Insert Here simulation Covariate
- }
- integG <-sum(eval(env$gFun,envir = env)*env$dt)
- env$s <- unique(c(env$s,t_k))[-length(env$s)]
- dumt <- env$t
- num <- length(env$Kern)
- integKer <- 0
- for(j in c(1:length(dumt))){
- env$t <- dumt[j]
- dumKernInt <- 0
- for(i in c(1:num)){
- lab.dx <- env$All.labKern at var.dx[i]
- dumKernInt <- dumKernInt+sum(eval(env$Kern,envir=env)*diff(eval(env[[lab.dx]])))
- }
- integKer <- integKer + dumKernInt
- }
- NewTerm <- 0
- if(env$Condcovariate){
- ## Insert Her
- }
- my.env$utrue<-u
- (log(env$U)+ integG + integKer+NewTerm)^2
- }
+ xinitCount <- numeric(length(object at Ppr@counting.var))
+ names(xinitCount) <- object at Ppr@counting.var
+ for(i in c(1:length(xinitCount))){
+ xinitCount[i] <- tail(env[[object at Ppr@counting.var[i]]],n = 1)
+ }
+ xinit <- c(xinit,xinitCount)
+ if(part[length(part)]-part[1]!=0){
+ DumVarCov <- simulate(yuima,
+ true.parameter = true.parameter,
+ increment.L = LevIncr,
+ sampling = setSampling(Terminal = (part[length(part)]-part[1]),
+ n = dim(LevIncr)[2]),
+ xinit=xinit[yuima at solve.variable])
+ for(i in c(1:length(object at Ppr@covariates))){
+ VariableDum <- DumVarCov at data@original.data[,yuima at solve.variable==object at Ppr@covariates[i]]
+ assign(object at Ppr@covariates[i], as.numeric(VariableDum), envir = env)
+ }
+ }else{
+ for(i in c(1:length(object at Ppr@covariates))){
+ VariableDum <- xinit[yuima at solve.variable==object at Ppr@covariates[i]]
+ assign(object at Ppr@covariates[i], as.numeric(VariableDum), envir = env)
+ }
+ }
+ #Insert Here simulation Covariate
+ }
+ integG <-sum(eval(env$gFun,envir = env)*env$dt)
+ env$s <- unique(c(env$s,t_k))[-length(env$s)]
+ dumt <- env$t
+ num <- length(env$Kern)
+ integKer <- 0
+ for(j in c(1:length(dumt))){
+ env$t <- dumt[j]
+ dumKernInt <- 0
+ for(i in c(1:num)){
+ lab.dx <- env$All.labKern at var.dx[i]
+ dumKernInt <- dumKernInt+sum(eval(env$Kern,envir=env)*diff(eval(env[[lab.dx]])))
+ }
+ integKer <- integKer + dumKernInt
+ }
+ NewTerm <- 0
+ if(env$Condcovariate){
+ ## Insert Her
+ }
+ my.env$utrue<-u
+ (log(env$U)+ integG + integKer+NewTerm)^2
+ }
u <- numeric(length = numbCountVar)
@@ -261,10 +315,10 @@
rownames(X_mat) <- object at model@solve.variable
dummyX <- simulate(object at model, true.parameter = true.parameter,
- increment.W = if(is.matrix(W[,1:pos])){W[,1:pos]}else{t(as.matrix(W[,1:pos]))},
- increment.L = if(is.matrix(dL[,1:pos])){dL[,1:pos]}else{t(as.matrix(dL[,1:pos]))},
- sampling = setSampling(Terminal = t_1,
- n = t_1/sampling at delta))
+ increment.W = if(is.matrix(W[,1:pos])){W[,1:pos]}else{t(as.matrix(W[,1:pos]))},
+ increment.L = if(is.matrix(dL[,1:pos])){dL[,1:pos]}else{t(as.matrix(dL[,1:pos]))},
+ sampling = setSampling(Terminal = t_1,
+ n = t_1/sampling at delta))
X_mat[,1:pos] <- t(dummyX at data@original.data)[,-1]
t_jump <- t_1
@@ -323,7 +377,7 @@
while(my.env$utrue<oldt_1){
assign("U",runif(1),envir = my.env)
optim((t_1+2*my.env$dt),Fun2,method = "Nelder-Mead",
- env=my.env)$par
+ env=my.env)$par
u[i] <- as.numeric(my.env$utrue)
}
}
@@ -343,47 +397,47 @@
#if(t_1!=sampling at Terminal){
- pos <- min(pos,dim(L)[2])
- JUMP[namesContVarJump, pos] <- L[namesContVarJump, pos]
- dL[object at Ppr@counting.var,c((pos0+1):pos)]<-JUMP[object at Ppr@counting.var,c((pos0+1):pos)]
- aa<-setSampling(Terminal = (t_1-my.env$old_u),
- n = length((pos0+1):pos))
- dummyX <- simulate(object at model, true.parameter = true.parameter,
- increment.W = if(is.matrix(W[,(pos0+1):pos])){W[,(pos0+1):pos]}else{t(as.matrix(W[,(pos0+1):pos]))},
- increment.L = if(is.matrix(dL[,(pos0+1):pos])){dL[,(pos0+1):pos]}else{t(as.matrix(dL[,(pos0+1):pos]))},
- sampling = aa,
- xinit=X_mat[,(pos0)])
- X_mat[,(pos0+1):pos] <- t(dummyX at data@original.data)[,-1]
- if(length(object at Kernel@variable.Integral at var.dx)==1){
- Comulat.dx <- apply(t(X_mat[object at Kernel@variable.Integral at var.dx,
- c((pos0+1):pos)]), 1, diff)
- }else{
- Comulat.dx <- apply(t(X_mat[object at Kernel@variable.Integral at var.dx,
- c((pos0+1):pos)]), 2, diff)
- }
- if(!is.matrix(Comulat.dx)){
- Comulat.dx <-t(as.matrix(Comulat.dx))
- }
+ pos <- min(pos,dim(L)[2])
+ JUMP[namesContVarJump, pos] <- L[namesContVarJump, pos]
+ dL[object at Ppr@counting.var,c((pos0+1):pos)]<-JUMP[object at Ppr@counting.var,c((pos0+1):pos)]
+ aa<-setSampling(Terminal = (t_1-my.env$old_u),
+ n = length((pos0+1):pos))
+ dummyX <- simulate(object at model, true.parameter = true.parameter,
+ increment.W = if(is.matrix(W[,(pos0+1):pos])){W[,(pos0+1):pos]}else{t(as.matrix(W[,(pos0+1):pos]))},
+ increment.L = if(is.matrix(dL[,(pos0+1):pos])){dL[,(pos0+1):pos]}else{t(as.matrix(dL[,(pos0+1):pos]))},
+ sampling = aa,
+ xinit=X_mat[,(pos0)])
+ X_mat[,(pos0+1):pos] <- t(dummyX at data@original.data)[,-1]
+ if(length(object at Kernel@variable.Integral at var.dx)==1){
+ Comulat.dx <- apply(t(X_mat[object at Kernel@variable.Integral at var.dx,
+ c((pos0+1):pos)]), 1, diff)
+ }else{
+ Comulat.dx <- apply(t(X_mat[object at Kernel@variable.Integral at var.dx,
+ c((pos0+1):pos)]), 2, diff)
+ }
+ if(!is.matrix(Comulat.dx)){
+ Comulat.dx <-t(as.matrix(Comulat.dx))
+ }
- Index <- matrix(c(1:prod(object at Kernel@Integrand at dimIntegrand)),
+ Index <- matrix(c(1:prod(object at Kernel@Integrand at dimIntegrand)),
nrow = object at Kernel@Integrand at dimIntegrand[1],
ncol = object at Kernel@Integrand at dimIntegrand[2])
- assign(object at Kernel@variable.Integral at var.time,
- sampling at grid[[1]][c((pos0+1):(pos))],
+ assign(object at Kernel@variable.Integral at var.time,
+ sampling at grid[[1]][c((pos0+1):(pos))],
envir = my.env)
- assign(object at gFun@param at time.var, t_1, envir = my.env)
- for(i in c(1:object at Kernel@Integrand at dimIntegrand[2])){
+ assign(object at gFun@param at time.var, t_1, envir = my.env)
+ for(i in c(1:object at Kernel@Integrand at dimIntegrand[2])){
- assign(object at Kernel@variable.Integral at var.dx[i],
+ assign(object at Kernel@variable.Integral at var.dx[i],
as.numeric(Comulat.dx[,i]),
envir = my.env)
- }
- pos0<-pos
- assign("pos0", pos, envir = my.env)
- assign("old_u",t_1, envir = my.env)
+ }
+ pos0<-pos
+ assign("pos0", pos, envir = my.env)
+ assign("old_u",t_1, envir = my.env)
- #}
- }
+ #}
+ }
assign("L",dL,envir = my.env)
}
X_mat[namesContVarJump,pos]<-X_mat[namesContVarJump,pos]
@@ -401,28 +455,28 @@
mynewincr <- if(is.matrix(res.dum$X_mat)){t(as.matrix(apply(cbind(0,res.dum$X_mat),1,diff)))}else{apply(cbind(0,res.dum$X_mat),1,diff)}
interResMod <- simulate(object = dummy.mod,
- true.parameter = true.parameter,
- sampling = sampling,
- increment.L = mynewincr)
+ true.parameter = true.parameter,
+ sampling = sampling,
+ increment.L = mynewincr)
- resGfun<-new("yuima.Output",
- Output = object at gFun,
- yuima=setYuima(model=dummy.mod,sampling = sampling))
+ resGfun<-new("yuima.Map",
+ Output = object at gFun,
+ yuima=setYuima(model=dummy.mod,sampling = sampling))
interResGfun <- simulate(object = resGfun,
- true.parameter = true.parameter,
- sampling = sampling,
- increment.L = mynewincr)
+ true.parameter = true.parameter,
+ sampling = sampling,
+ increment.L = mynewincr)
dummyObject <- object at Kernel
dummyObject at variable.Integral@out.var <-object at Ppr@additional.info
resInt <- new("yuima.Integral",
- Integral = dummyObject,
- yuima = setYuima(model=dummy.mod,sampling = sampling))
+ Integral = dummyObject,
+ yuima = setYuima(model=dummy.mod,sampling = sampling))
interResInt <- simulate(object = resInt,
- true.parameter = true.parameter,
- sampling = sampling,
- increment.L = mynewincr)
+ true.parameter = true.parameter,
+ sampling = sampling,
+ increment.L = mynewincr)
DataIntensity <- interResGfun at data@original.data + interResInt at data@original.data
InterMDia<-zoo(interResMod at data@original.data, order.by = index(DataIntensity))
Alldata <-merge(InterMDia,DataIntensity)
@@ -437,11 +491,9 @@
# }
object at data<-setData(Alldata)
return(object)
- }
+}
-
-
# simOzaki.aux<-function(gFun,a,cCoeff, Time, numJump){
# t_k<-0
# N<-0
More information about the Yuima-commits
mailing list