[Yuima-commits] r573 - pkg/yuima/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Jan 22 09:59:49 CET 2017
Author: lorenzo
Date: 2017-01-22 09:59:49 +0100 (Sun, 22 Jan 2017)
New Revision: 573
Modified:
pkg/yuima/R/AuxMethodforPPR.R
Log:
QMLE for univariate PPR
Modified: pkg/yuima/R/AuxMethodforPPR.R
===================================================================
--- pkg/yuima/R/AuxMethodforPPR.R 2017-01-22 06:10:17 UTC (rev 572)
+++ pkg/yuima/R/AuxMethodforPPR.R 2017-01-22 08:59:49 UTC (rev 573)
@@ -2,141 +2,209 @@
## Regression Model
is.Ppr <- function(yuimaPpr){is(yuimaPpr,"yuima.Ppr")}
+Internal.LogLikPPR <- function(param,my.envd1=NULL,
+ my.envd2=NULL,my.envd3=NULL){
+ param<-unlist(param)
+
+ IntLambda<-InternalConstractionIntensity(param,my.envd1,
+ my.envd2,my.envd3)
+ Index<-my.envd3$gridTime
+ Integr1 <- -sum(IntLambda[-length(IntLambda)]*diff(Index))
+ cond2 <- diff(as.numeric(my.envd3$YUIMA.PPR at data@original.data))
+ Integr2<- sum(log(IntLambda[-1][cond2>0]))
+ logLik <- Integr1+Integr2
+ cat("\n ",logLik)
+ return(-logLik)
+}
+
quasiLogLik.Ppr <- function(yuimaPpr, parLambda=list(), method=method, fixed = list(),
- lower, upper, call, ...){
- PprData<-yuimaPpr at data
- Time <- index(yuimaPpr at data@zoo.data[[1]])
- envPpr <- list()
- dY <- paste0("d",yuimaPpr at Ppr@var.dx)
- for(i in (c(1:(length(Time)-1)))){
- envPpr[[i]]<-new.env()
- assign(yuimaPpr at gFun@param at time.var,rep(Time[i+1],i),envir=envPpr[[i]])
- assign(yuimaPpr at Ppr@var.dt,Time[1:i],envir=envPpr[[i]])
- if(length(yuimaPpr at Ppr@covariates)>0){
- for(j in c(1:length(yuimaPpr at Ppr@covariates))){
- cond<-colnames(yuimaPpr at data@original.data)%in%yuimaPpr at Ppr@covariates[[j]]
- assign(yuimaPpr at Ppr@covariates[[j]],
- as.numeric(yuimaPpr at data@original.data[1:(i+1),cond]),
- envir=envPpr[[i]])
+ lower, upper, call, ...){
+
+ yuimaPpr->yuimaPPr
+ parLambda->param
+ gfun<-yuimaPPr at gFun@formula
+
+ dimIntegr <- length(yuimaPPr at Kernel@Integrand at IntegrandList)
+ Integrand2 <- character(length=dimIntegr)
+ for(i in c(1:dimIntegr)){
+ Integrand1 <- as.character(yuimaPPr at Kernel@Integrand at IntegrandList[[i]])
+ timeCond <- paste0(" * (",yuimaPPr at Kernel@variable.Integral at var.time," < ",yuimaPPr at Kernel@variable.Integral at upper.var,")")
+ Integrand2[i] <-paste0(Integrand1,timeCond)
+ }
+
+ Integrand2<- matrix(Integrand2,yuimaPPr at Kernel@Integrand at dimIntegrand[1],yuimaPPr at Kernel@Integrand at dimIntegrand[2])
+
+
+ for(j in c(1:yuimaPPr at Kernel@Integrand at dimIntegrand[2])){
+ Integrand2[,j]<-paste0(Integrand2[,j]," * d",yuimaPPr at Kernel@variable.Integral at var.dx[j])
+ }
+ colnames(Integrand2) <- paste0("d",yuimaPPr at Kernel@variable.Integral at var.dx)
+ NamesIntegrandExpr <- as.character(matrix(colnames(Integrand2), dim(Integrand2)[1],dim(Integrand2)[2], byrow = TRUE))
+ Integrand2expr<- parse(text=Integrand2)
+
+ gridTime <- time(yuimaPPr at data@original.data)
+
+ yuimaPPr at Kernel@variable.Integral at var.dx
+ if(any(yuimaPPr at Kernel@variable.Integral at var.dx %in% yuimaPPr at model@solve.variable)){
+ my.envd1<-new.env()
+ ExistdN<-TRUE
+ }else{
+ ExistdN<-FALSE
+ }
+ Univariate<-FALSE
+ if(length(yuimaPPr at Ppr@counting.var)==1){
+ Univariate<-TRUE
+ }
+ if(any(!(yuimaPPr at Kernel@variable.Integral at var.dx %in% yuimaPPr at model@solve.variable))){
+ my.envd2<-new.env()
+ ExistdX<-TRUE
+ }else{
+ my.envd2<-new.env()
+ ExistdX<-FALSE
+ }
+
+ my.envd3 <- new.env()
+ namesparam<-names(param)
+ if(!(all(namesparam %in% yuimaPPr at Ppr@allparamPpr) && length(namesparam)==length(yuimaPPr at Ppr@allparamPpr))){
+ return(NULL)
+ }
+
+ # construction my.envd1
+ if(ExistdN){
+
+ #CountingVariable
+ for(i in c(1:length(yuimaPPr at Ppr@counting.var))){
+ cond <- yuimaPPr at Ppr@counting.var[i] %in% yuimaPPr at model@solve.variable
+ dummyData <-unique(yuimaPPr at data@original.data[,cond])[-1]
+ assign(yuimaPPr at Ppr@counting.var[i], dummyData,envir=my.envd1)
+ }
+ # Names expression
+ assign("NamesIntgra", NamesIntegrandExpr, envir=my.envd1)
+ #dN
+ namedX <-NULL
+ for(i in c(1:length(yuimaPPr at Kernel@variable.Integral at var.dx))){
+ if(yuimaPPr at Kernel@variable.Integral at var.dx[i] %in% yuimaPPr at Ppr@counting.var){
+ cond <- yuimaPPr at model@solve.variable %in% yuimaPPr at Kernel@variable.Integral at var.dx[i]
+ namedX<-c(namedX,paste0("d",yuimaPPr at Kernel@variable.Integral at var.dx[i]))
+ dummyData <- diff(as.numeric(yuimaPPr at data@original.data[,cond]))# We consider only Jump
+ dummyJumpTime <- gridTime[-1][dummyData>0]
+ dummyData2 <- diff(unique(cumsum(dummyData)))
+ dummyData3 <- zoo(dummyData2,order.by = dummyJumpTime)
+ assign(paste0("d",yuimaPPr at Kernel@variable.Integral at var.dx[i]), dummyData3 ,envir=my.envd1)
}
}
- for(j in c(1:length(yuimaPpr at Ppr@counting.var))){
- #cond<-colnames(yuimaPpr at data@original.data)%in%yuimaPpr at Ppr@counting.var[[j]]
- cond<-yuimaPpr at model@solve.variable%in%yuimaPpr at Ppr@counting.var[[j]]
- assign(yuimaPpr at Ppr@counting.var[[j]],
- as.numeric(yuimaPpr at data@original.data[1:(i+1),cond]),
- envir=envPpr[[i]])
+ assign("namedX",namedX, envir = my.envd1)
+ assign("var.time",yuimaPPr at Kernel@variable.Integral at var.time,envir=my.envd1)
+ assign("t.time",yuimaPPr at Kernel@variable.Integral at upper.var,envir=my.envd1)
+
+ # Covariates
+ if(length(yuimaPPr at Ppr@covariates)>1){
+ # Covariates should be identified at jump time
+ return(NULL)
}
- for(j in c(1:length(yuimaPpr at Ppr@var.dx))){
- #cond<-c(colnames(yuimaPpr at data@original.data),yuimaPpr at Ppr@var.dt)%in%c(yuimaPpr at Ppr@var.dx,yuimaPpr at Ppr@var.dt)[[j]]
- cond<-c(yuimaPpr at model@solve.variable,yuimaPpr at Ppr@var.dt)%in%c(yuimaPpr at Ppr@var.dx,yuimaPpr at Ppr@var.dt)[[j]]
- if(any(cond[-length(cond)])){
- assign(paste0("d",yuimaPpr at Ppr@var.dx[[j]]),
- diff(as.numeric(yuimaPpr at data@original.data[1:(i+1),cond[-length(cond)]])),
- envir=envPpr[[i]])
- }
- if(tail(cond,n=1L)){
- assign(paste0("d",yuimaPpr at Ppr@var.dx[[j]]),
- as.numeric(diff(Time[1:(i+1)])),
- envir=envPpr[[i]])
- }
- }
}
- IntKernExpr<- function(Kern,dy){
- dum<-paste(Kern,dy,sep="*")
- dum<-paste0(dum, collapse = " + ")
- dum <- paste0("sum( ", dum, " )")
- return(parse(text = dum))
- }
+ # end coonstruction my.envd1
- IntegKern <- lapply(yuimaPpr at Kernel@Integrand at IntegrandList,IntKernExpr,dY)
- Integrator <- t(as.matrix(eval(parse(text=dY[1]),envir=envPpr[[length(envPpr)]])))
- if(length(dY)>1){
- for(i in c(2:length(dY))){
- Integrator <- rbind(Integrator,
- t(as.matrix(eval(parse(text=dY[1]),envir=envPpr[[length(envPpr)]]))))
+ # construction my.envd2
+ if(ExistdX){
+ #Covariate
+
+ #CountingVariable
+ for(i in c(1:length(yuimaPPr at Ppr@counting.var))){
+ cond <- yuimaPPr at Ppr@counting.var[i] %in% yuimaPPr at model@solve.variable
+ dummyData <-yuimaPPr at data@original.data[,cond]
+ assign(yuimaPPr at Ppr@counting.var[i], dummyData,envir=my.envd1)
}
+
+
+ }else{
+ assign("KerneldX",NULL,envir=my.envd2)
}
- assign("Integrator",Integrator,envir=envPpr[[length(envPpr)]])
- assign("Nlamb",length(yuimaPpr at Ppr@counting.var),envir=envPpr[[length(envPpr)]])
- cond <- yuimaPpr at model@solve.variable %in% yuimaPpr at Ppr@counting.var
- #assign("CountVar",yuimaPpr at data@original.data[,cond],envir=envPpr[[length(envPpr)]])
- DummyCountVar <- yuimaPpr at data@original.data[,cond]
- condDummyCountVar <- t(apply(as.matrix(DummyCountVar),FUN = "diff",MARGIN = 2))!=0
- assign("CountVar", condDummyCountVar, envir=envPpr[[length(envPpr)]])
+ # end construction my.envd2
+
+ # construction my.envd3
+
+ #Covariate
+
+ #CountingVariable
+ for(i in c(1:length(yuimaPPr at Ppr@counting.var))){
+ cond <- yuimaPPr at Ppr@counting.var[i] %in% yuimaPPr at model@solve.variable
+ dummyData <-yuimaPPr at data@original.data[,cond]
+ assign(yuimaPPr at Ppr@counting.var[i], dummyData,envir=my.envd3)
+ }
+ #time
+ assign(yuimaPPr at model@time.variable, gridTime, my.envd3)
+
+ #Model
+ assign("YUIMA.PPR",yuimaPPr,envir=my.envd3)
+ assign("namesparam",namesparam,envir=my.envd3)
+ assign("gfun",gfun,envir=my.envd3)
+ assign("Integrand2",Integrand2,envir=my.envd3)
+ assign("Integrand2expr",Integrand2expr,envir=my.envd3)
+
+ assign("gridTime",gridTime,envir=my.envd3)
+ assign("Univariate",Univariate,envir=my.envd3)
+ assign("ExistdN",ExistdN,envir=my.envd3)
+ assign("ExistdX",ExistdX,envir=my.envd3)
out<-NULL
- param1<-unlist(parLambda)
- my.env <- envPpr
+
if(length(lower)==0 && length(upper)>0 && length(fixed)==0){
- out <- optim(par=param1, fn=aux.lambdaFromData,
- method = method, upper=upper, envPpr = my.env, gFun=yuimaPpr at gFun,
- Kern =IntegKern, intensityParm = yuimaPpr at Ppr@allparamPpr,
- logLikelihood=TRUE, ...)
+ out <- optim(par=param, fn=Internal.LogLikPPR,
+ my.envd1=my.envd1,my.envd2=my.envd2,my.envd3=my.envd3,
+ method = method, upper=upper, ...)
}
if(length(lower)==0 && length(upper)==0 && length(fixed)>0){
- out <- optim(par = param1, fn = aux.lambdaFromData,
- method = method, fixed = fixed, envPpr = my.env,
- gFun = yuimaPpr at gFun, Kern = IntegKern,
- intensityParm = yuimaPpr at Ppr@allparamPpr,
- logLikelihood = TRUE, ...)
+ out <- optim(par=param, fn=Internal.LogLikPPR,
+ my.envd1=my.envd1,my.envd2=my.envd2,my.envd3=my.envd3,
+ method = method, fixed = fixed, ...)
}
if(length(lower)>0 && length(upper)==0 && length(fixed)==0){
- out <- optim(par = param1, fn = aux.lambdaFromData,
- method = method, lower=lower, envPpr = my.env,
- gFun = yuimaPpr at gFun, Kern = IntegKern,
- intensityParm = yuimaPpr at Ppr@allparamPpr,
- logLikelihood = TRUE, ...)
+ out <- optim(par = param, fn=Internal.LogLikPPR,
+ my.envd1=my.envd1,my.envd2=my.envd2,my.envd3=my.envd3,
+ method = method, lower=lower, ...)
}
if(length(lower)>0 && length(upper)>0 && length(fixed)==0){
- out <- optim(par = param1, fn = aux.lambdaFromData,
+ out <- optim(par=param, fn=Internal.LogLikPPR,
+ my.envd1=my.envd1,my.envd2=my.envd2,my.envd3=my.envd3,
method = method, upper = upper,
- lower=lower, envPpr = my.env,
- gFun = yuimaPpr at gFun, Kern = IntegKern,
- intensityParm = yuimaPpr at Ppr@allparamPpr,
- logLikelihood = TRUE, ...)
+ lower=lower, ...)
}
if(length(lower)==0 && length(upper)>0 && length(fixed)>0){
- out <- optim(par = param1, fn = aux.lambdaFromData,
+ out <- optim(par=param, fn=Internal.LogLikPPR,
+ my.envd1=my.envd1,my.envd2=my.envd2,my.envd3=my.envd3,
method = method, upper = upper,
- fixed = fixed, envPpr = my.env,
- gFun = yuimaPpr at gFun, Kern = IntegKern,
- intensityParm = yuimaPpr at Ppr@allparamPpr,
- logLikelihood = TRUE, ...)
+ fixed = fixed, ...)
}
if(length(lower)>0 && length(upper)==0 && length(fixed)>0){
- out <- optim(par=param1, fn=aux.lambdaFromData,
+ out <- optim(par=param, fn=Internal.LogLikPPR,
+ my.envd1=my.envd1,my.envd2=my.envd2,my.envd3=my.envd3,
method = method, lower = lower,
- fixed = fixed, envPpr = my.env,
- gFun = yuimaPpr at gFun, Kern = IntegKern,
- intensityParm = yuimaPpr at Ppr@allparamPpr,
- logLikelihood = TRUE, ...)
+ fixed = fixed, ...)
}
if(length(lower)>0 && length(upper)>0 && length(fixed)>0){
- out <- optim(par = param1, fn = aux.lambdaFromData,
- method = method, lower = lower, fixed = fixed, upper = upper,
- envPpr = my.env, gFun = yuimaPpr at gFun, Kern = IntegKern,
- intensityParm = yuimaPpr at Ppr@allparamPpr, logLikelihood = TRUE, ...)
+ out <- optim(par=param, fn=Internal.LogLikPPR,
+ my.envd1=my.envd1,my.envd2=my.envd2,my.envd3=my.envd3,
+ method = method, lower = lower, fixed = fixed, upper = upper, ...)
}
if(is.null(out)){
- out <- optim(par = param1, fn = aux.lambdaFromData,
- method = method, envPpr = my.env, gFun = yuimaPpr at gFun,
- Kern = IntegKern, intensityParm = yuimaPpr at Ppr@allparamPpr,
- logLikelihood = TRUE, ...)
+ out <- optim(par=param, fn=Internal.LogLikPPR,
+ my.envd1=my.envd1,my.envd2=my.envd2,my.envd3=my.envd3,
+ method = method, ...)
}
@@ -144,6 +212,7 @@
}
+
lambdaFromData <- function(yuimaPpr, PprData=NULL, parLambda=list()){
if(is.null(PprData)){
PprData<-yuimaPpr at data
More information about the Yuima-commits
mailing list