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

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Thu Mar 21 20:59:51 CET 2019


Author: lorenzo
Date: 2019-03-21 20:59:50 +0100 (Thu, 21 Mar 2019)
New Revision: 690

Modified:
   pkg/yuima/R/AuxMethodforPPR.R
   pkg/yuima/R/simulateMultiProcess.R
Log:


Modified: pkg/yuima/R/AuxMethodforPPR.R
===================================================================
--- pkg/yuima/R/AuxMethodforPPR.R	2019-03-04 19:20:51 UTC (rev 689)
+++ pkg/yuima/R/AuxMethodforPPR.R	2019-03-21 19:59:50 UTC (rev 690)
@@ -78,7 +78,7 @@
 
   yuimaPPR->yuimaPPR
   parLambda->param
-  gfun<-yuimaPPR at gFun@formula
+ # gfun<-yuimaPPR at gFun@formula
 
   gfun<-yuimaPPR at gFun@formula
 
@@ -112,7 +112,7 @@
 
   gridTime <- time(yuimaPPR at data@original.data)
 
-  yuimaPPR at Kernel@variable.Integral at var.dx
+ # 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
@@ -134,46 +134,55 @@
   my.envd3 <- new.env()
   namesparam<-names(param)
   resCov<-NULL
+  NoFeedBackIntensity <- TRUE
   if(!(all(namesparam %in% yuimaPPR at PPR@allparamPPR) && length(namesparam)==length(yuimaPPR at PPR@allparamPPR))){
     
     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]])))
+      
+      if(!all(yuimaPPR at model@state.variable %in% yuimaPPR at model@solve.variable)|yuimaPPR at PPR@RegressWithCount){
+        NoFeedBackIntensity <- FALSE
       }
-      
-      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
+      if(NoFeedBackIntensity){  
+        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)){
-          dummyJump0 <- c(dummyJump0,
-                          as.character(unlist(yuimaPPR at model@jump.coeff[[j]][1:dimJump])))
+          dummydiff0<-c(dummydiff0,
+                        as.character(unlist(yuimaPPR at model@diffusion[[j]])))
         }
-        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)
+        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)
     }      
   }
 
@@ -439,9 +448,76 @@
                   method = method, nobs=as.integer(N.jump), model=my.envd3$YUIMA.PPR)
   if(!is.null(resCov)){
    return(list(PPR=final_res,Covariates=resCov))
-  }else{
+  }
+  if(!NoFeedBackIntensity){
+    if(all(yuimaPPR at model@state.variable %in% yuimaPPR at model@solve.variable)){
+      myMod <- yuimaPPR at model
+      myYuima <- setYuima(data = yuimaPPR at data, 
+        model = yuimaPPR at model, sampling = yuimaPPR at sampling)
+      resCov <- qmleLevy(yuima = myYuima,
+        start=param[myMod at parameter@all],upper=upper[myMod at parameter@all],
+        lower=lower[myMod at parameter@all])
+    }else{
+      IntensityData <- Intensity.PPR(final_res at model,
+                                     param=coef(final_res))
+      mylambda <- IntensityData at original.data
+      OrigData <- yuimaPPR at data@original.data
+      NewData0 <- cbind(OrigData,mylambda)
+      colnames(NewData0) <- yuimaPPR at model@state.variable
+      NewData<-setData(zoo(NewData0,
+                           order.by = index(IntensityData at zoo.data[[1]])))
+      
+      lengthVar <- length(yuimaPPR at model@state.variable)
+      lengthOrigVar <- length(yuimaPPR at model@solve.variable)
+      DummyDrift <- as.character(rep(0,lengthVar))
+      DummyDrift[1:lengthOrigVar] <- as.character(yuimaPPR at model@drift) 
+      
+      dummydiff0<- NULL
+      for(j in c(1:lengthOrigVar)){
+        dummydiff0<-c(dummydiff0,
+                      as.character(unlist(yuimaPPR at model@diffusion[[j]])))
+      }
+      
+      dummydiff <- matrix(dummydiff0, nrow = lengthOrigVar, 
+                          ncol = length(dummydiff0)/lengthOrigVar)
+      
+      if(length(yuimaPPR at model@jump.variable)!=0){
+        dummydiff <- rbind(dummydiff,matrix("0",
+                                            nrow = lengthVar-lengthOrigVar,dim(dummydiff)[2]))
+      }
+      
+      dummyJump0 <- NULL
+      for(j in c(1:lengthOrigVar)){
+        dummyJump0 <- c(dummyJump0,
+                        as.character(unlist(yuimaPPR at model@jump.coeff[[j]][])))
+      }
+      dummyJump <- matrix(dummyJump0, nrow=lengthOrigVar,ncol=length(dummyJump0)/lengthOrigVar, byrow = T)
+      
+      if(length(yuimaPPR at model@jump.variable)!=0){
+        dummyJump1<- matrix(as.character(diag(lengthVar)),lengthVar,lengthVar) 
+        dummyJump1[1:lengthOrigVar,1:lengthOrigVar] <- dummyJump
+      }  
+      
+      # aaa<-setModel(drift="1",diffusion = "1")
+      # aaa at jump.variable
+      # yuimaPPR at model@jump.variable
+      meas.type <- rep("code",lengthVar)
+      myMod <- setModel(drift = DummyDrift, diffusion = dummydiff,
+                        jump.coeff = dummyJump1, jump.variable = yuimaPPR at model@jump.variable,
+                        measure = list(df=yuimaPPR at model@measure$df),
+                        measure.type = meas.type,
+                        solve.variable = yuimaPPR at model@state.variable,
+                        state.variable = yuimaPPR at model@state.variable)
+      myYuima <- setYuima(data = NewData, model = myMod)
+      resCov <- qmleLevy(yuima = myYuima,
+                      start=param[myMod at parameter@all],upper=upper[myMod at parameter@all],
+                      lower=lower[myMod at parameter@all])
+    }
+    return(list(PPR=final_res,Covariates=resCov))
+  } 
+  
    return(final_res)
-  }
+  
 }
 
 

Modified: pkg/yuima/R/simulateMultiProcess.R
===================================================================
--- pkg/yuima/R/simulateMultiProcess.R	2019-03-04 19:20:51 UTC (rev 689)
+++ pkg/yuima/R/simulateMultiProcess.R	2019-03-21 19:59:50 UTC (rev 690)
@@ -335,72 +335,73 @@
           #yuima.stop("code multivariate Levy will be implemented as soon as possible")
         }
      }else{
-       intens <- as.character(sdeModel at measure$intensity)
-       dens <- as.character(sdeModel at measure$df$expr)
-       # If we consider independence between CP and the Other Levy
-       # we have:
-       numbLev <- length(sdeModel at measure.type)
-       posCPindex <- c(1:numbLev)[sdeModel at measure.type%in%"CP"]
-       CPmeasureComp <- paste0(dens,"[,c(",toString(posCPindex),")]")
-       intens <- as.character(sdeModel at measure$intensity)
-       dumCP <- setPoisson(intensity = intens, df = CPmeasureComp,
-         dimension = length(posCPindex))
-
-       # Simulation CP part
-       dummSamp <- yuima at sampling
-       samp <- setSampling(Initial = dummSamp at Initial,
-                           Terminal = unique(dummSamp at Terminal),
-                           n = unique(dummSamp at n))
-       trajCP <- simulate(object = dumCP, sampling = samp,
-          true.parameter = true.parameter)
-
-
-       dimJumpCoeff <- length(yuima at model@measure.type)
-       dumjumpCoeff <- matrix(as.character(diag(rep(1,dimJumpCoeff))),dimJumpCoeff,dimJumpCoeff)
-       Dumsolve.variable <- paste0("MyLevyDum",c(1:dimJumpCoeff))
-       dummy.measure.code <- as.character(sdeModel at measure$df$expr)
-       LevyMod <- setMultiModel(drift=rep("0",dimJumpCoeff),
-                                diffusion = NULL,
-                                jump.coeff = dumjumpCoeff,
-                                df = dummy.measure.code,
-                                measure.type = "code",
-                                solve.variable = Dumsolve.variable)
-       yuimaLevy <- setYuima(model=LevyMod, sampling = samp)
-       yuimaLevy at model@dimension <- dimJumpCoeff
-
-       trajcode<- simCode(xinit=rep("0",length=dimJumpCoeff),
-          yuima = yuimaLevy, env=yuimaEnv)
-
-       countCP <- 0
-       countcode <- 0
-       if(yuima at model@measure.type[1]=="CP"){
-          Incr.levy <- as.matrix(as.numeric(diff(trajCP at data@zoo.data[[1]])))
-          countcode <- countcode+1
-       }else{
-         if(yuima at model@measure.type[1]=="code"){
-            Incr.levy <- as.matrix(as.numeric(diff(trajcode at data@zoo.data[[1]])))
-            countCP <- countCP+1
+       if(any(sdeModel at measure.type=="CP")){
+         intens <- as.character(sdeModel at measure$intensity)
+         dens <- as.character(sdeModel at measure$df$expr)
+         # If we consider independence between CP and the Other Levy
+         # we have:
+         numbLev <- length(sdeModel at measure.type)
+         posCPindex <- c(1:numbLev)[sdeModel at measure.type%in%"CP"]
+         CPmeasureComp <- paste0(dens,"[,c(",toString(posCPindex),")]")
+         intens <- as.character(sdeModel at measure$intensity)
+         dumCP <- setPoisson(intensity = intens, df = CPmeasureComp,
+           dimension = length(posCPindex))
+  
+         # Simulation CP part
+         dummSamp <- yuima at sampling
+         samp <- setSampling(Initial = dummSamp at Initial,
+                             Terminal = unique(dummSamp at Terminal),
+                             n = unique(dummSamp at n))
+         trajCP <- simulate(object = dumCP, sampling = samp,
+            true.parameter = true.parameter)
+  
+  
+         dimJumpCoeff <- length(yuima at model@measure.type)
+         dumjumpCoeff <- matrix(as.character(diag(rep(1,dimJumpCoeff))),dimJumpCoeff,dimJumpCoeff)
+         Dumsolve.variable <- paste0("MyLevyDum",c(1:dimJumpCoeff))
+         dummy.measure.code <- as.character(sdeModel at measure$df$expr)
+         LevyMod <- setMultiModel(drift=rep("0",dimJumpCoeff),
+                                  diffusion = NULL,
+                                  jump.coeff = dumjumpCoeff,
+                                  df = dummy.measure.code,
+                                  measure.type = "code",
+                                  solve.variable = Dumsolve.variable)
+         yuimaLevy <- setYuima(model=LevyMod, sampling = samp)
+         yuimaLevy at model@dimension <- dimJumpCoeff
+  
+         trajcode<- simCode(xinit=rep("0",length=dimJumpCoeff),
+            yuima = yuimaLevy, env=yuimaEnv)
+  
+         countCP <- 0
+         countcode <- 0
+         if(yuima at model@measure.type[1]=="CP"){
+            Incr.levy <- as.matrix(as.numeric(diff(trajCP at data@zoo.data[[1]])))
+            countcode <- countcode+1
+         }else{
+           if(yuima at model@measure.type[1]=="code"){
+              Incr.levy <- as.matrix(as.numeric(diff(trajcode at data@zoo.data[[1]])))
+              countCP <- countCP+1
+           }
          }
-       }
-
-       if(length(yuima at model@measure.type)>1){
-         for(i in c(2:length(yuima at model@measure.type))){
-           if(yuima at model@measure.type[i]=="CP"){
-             Incr.levy<-cbind(Incr.levy,as.numeric(diff(trajCP at data@zoo.data[[(i-countCP)]])))
-             countcode <- countcode+1
-           }else{
-            if(yuima at model@measure.type[i]=="code"){
-              Incr.levy <- cbind(Incr.levy,as.numeric(diff(trajcode at data@zoo.data[[i]])))
-              countCP <- countCP+1
+  
+         if(length(yuima at model@measure.type)>1){
+           for(i in c(2:length(yuima at model@measure.type))){
+             if(yuima at model@measure.type[i]=="CP"){
+               Incr.levy<-cbind(Incr.levy,as.numeric(diff(trajCP at data@zoo.data[[(i-countCP)]])))
+               countcode <- countcode+1
+             }else{
+              if(yuima at model@measure.type[i]=="code"){
+                Incr.levy <- cbind(Incr.levy,as.numeric(diff(trajcode at data@zoo.data[[i]])))
+                countCP <- countCP+1
+              }
+             }
             }
-           }
-          }
+         }
+  
+  
+          # yuima.stop("Levy with CP and/or code")
        }
-
-
-        # yuima.stop("Levy with CP and/or code")
      }
-
    }
    if(!is.null(increment.L))
      Incr.levy<-t(increment.L)



More information about the Yuima-commits mailing list