[Yuima-commits] r483 - in pkg/yuimaGUI: . inst/yuimaGUI inst/yuimaGUI/www

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Sun Oct 16 21:05:53 CEST 2016


Author: phoenix844
Date: 2016-10-16 21:05:53 +0200 (Sun, 16 Oct 2016)
New Revision: 483

Modified:
   pkg/yuimaGUI/DESCRIPTION
   pkg/yuimaGUI/inst/yuimaGUI/global.R
   pkg/yuimaGUI/inst/yuimaGUI/server.R
   pkg/yuimaGUI/inst/yuimaGUI/ui.R
   pkg/yuimaGUI/inst/yuimaGUI/www/custom.css
Log:
Added yuima cpoint estimation. Some troubles with carma and cogarch, don't use them

Modified: pkg/yuimaGUI/DESCRIPTION
===================================================================
--- pkg/yuimaGUI/DESCRIPTION	2016-10-15 15:14:36 UTC (rev 482)
+++ pkg/yuimaGUI/DESCRIPTION	2016-10-16 19:05:53 UTC (rev 483)
@@ -1,7 +1,7 @@
 Package: yuimaGUI
 Type: Package
 Title: A Graphical User Interface for the Yuima Package
-Version: 0.7.4
+Version: 0.7.5
 Author: YUIMA Project Team
 Maintainer: Emanuele Guidotti <emanuele.guidotti at studenti.unimi.it>
 Description: Provides a graphical user interface for the yuima package.

Modified: pkg/yuimaGUI/inst/yuimaGUI/global.R
===================================================================
--- pkg/yuimaGUI/inst/yuimaGUI/global.R	2016-10-15 15:14:36 UTC (rev 482)
+++ pkg/yuimaGUI/inst/yuimaGUI/global.R	2016-10-16 19:05:53 UTC (rev 483)
@@ -8,12 +8,13 @@
 require(shinyBS)
 #require(corrplot)
 
+options(warn=-1)
 
 if(!exists("yuimaGUItable"))
   yuimaGUItable <<- reactiveValues(series=data.frame(),  model=data.frame(), simulation=data.frame(), hedging=data.frame())
 
 if(!exists("yuimaGUIdata"))
-  yuimaGUIdata <<- reactiveValues(series=list(), cp=list(), model=list(), simulation=list(), hedging = list())
+  yuimaGUIdata <<- reactiveValues(series=list(), cp=list(), cpYuima=list(), model=list(), simulation=list(), hedging = list())
 
 if(!exists("estimateSettings"))
   estimateSettings <<- list()
@@ -124,9 +125,9 @@
 
 observe({
   differ <- names(yuimaGUIdata$cp)[!(names(yuimaGUIdata$cp) %in% names(yuimaGUIdata$series))]
-  if (length(differ)!=0)
-    for (i in differ)
-      yuimaGUIdata$cp[[i]] <<- NULL
+  if (length(differ)!=0) for (i in differ) yuimaGUIdata$cp[[i]] <<- NULL
+  differ <- names(yuimaGUIdata$cpYuima)[!(names(yuimaGUIdata$cpYuima) %in% names(yuimaGUIdata$series))]
+  if (length(differ)!=0) for (i in differ) yuimaGUIdata$cpYuima[[i]] <<- NULL
 })
 
 
@@ -194,7 +195,7 @@
                     "Compound Poisson" = "Power Low Intensity",
                     "Compound Poisson" = "Exponentially Decaying Intensity",
                     "Compound Poisson" = "Periodic Intensity",
-                    "CARMA" = "Carma Noise: Compound Poisson",
+                    "CARMA" = "Carma(p,q)",
                     "COGARCH" = "Cogarch(p,q)"
                     )
 
@@ -225,29 +226,29 @@
     names(startmax) <- par
     startmin["a0"] <- ifelse(is.na(lower),NA,0)
     startmax["a0"] <- ifelse(is.na(upper),NA,1000)
-    if (!is.na(jumps)){
-      boundsJump <- jumpBounds(jumps = jumps, lower = lower, upper = upper)
-      for (i in par[par %in% names(boundsJump$lower)]){
-        startmin[[i]] <- boundsJump$lower[[i]]
-        startmax[[i]] <- boundsJump$upper[[i]]
-      }
-    }
+#     if (!is.na(jumps)){
+#       boundsJump <- jumpBounds(jumps = jumps, lower = lower, upper = upper)
+#       for (i in par[par %in% names(boundsJump$lower)]){
+#         startmin[[i]] <- boundsJump$lower[[i]]
+#         startmax[[i]] <- boundsJump$upper[[i]]
+#       }
+#     }
     return(list(lower=as.list(startmin), upper=as.list(startmax)))
   }
   if (name %in% defaultModels[names(defaultModels) == "CARMA"]){
     par <- setModelByName(name = name, jumps = jumps,  AR_C = AR_C, MA_C = MA_C)@parameter
-    par <- unique(c(par at drift, par at jump, par at measure, par at xinit[1]))
+    par <- unique(c(par at drift, par at xinit))
     startmin <- rep(lower, length(par))
     startmax <- rep(upper, length(par))
     names(startmin) <- par
     names(startmax) <- par
-    if (!is.na(jumps)){
-      boundsJump <- jumpBounds(jumps = jumps, lower = lower, upper = upper)
-      for (i in par[par %in% names(boundsJump$lower)]){
-        startmin[[i]] <- boundsJump$lower[[i]]
-        startmax[[i]] <- boundsJump$upper[[i]]
-      }
-    }
+#     if (!is.na(jumps)){
+#       boundsJump <- jumpBounds(jumps = jumps, lower = lower, upper = upper)
+#       for (i in par[par %in% names(boundsJump$lower)]){
+#         startmin[[i]] <- boundsJump$lower[[i]]
+#         startmax[[i]] <- boundsJump$upper[[i]]
+#       }
+#     }
     return(list(lower=as.list(startmin), upper=as.list(startmax)))
   }
   if (name == "Brownian Motion" | name == "Bm")
@@ -327,7 +328,7 @@
 }
 
 
-setModelByName <- function(name, jumps, AR_C = NA, MA_C = NA, XinExpr = FALSE){
+setModelByName <- function(name, jumps = NA, AR_C = NA, MA_C = NA, XinExpr = FALSE){
   if (name %in% names(isolate({usr_models$model}))){
     if (isolate({usr_models$model[[name]]$class=="Diffusion process"}))
       return(isolate({usr_models$model[[name]]$object}))
@@ -357,15 +358,8 @@
   if (name == "Linear Intensity") return(yuima::setPoisson(intensity="alpha+beta*t", df=setJumps(jumps = jumps), solve.variable = "x"))
   if (name == "Exponentially Decaying Intensity") return(yuima::setPoisson(intensity="alpha*exp(-beta*t)", df=setJumps(jumps = jumps), solve.variable = "x"))
   if (name == "Periodic Intensity") return(yuima::setPoisson(intensity="a/2*(1+cos(omega*t+phi))+b", df=setJumps(jumps = jumps), solve.variable = "x"))
-  if (name == "Cogarch(p,q)") {
-    return(yuima::setCogarch(p = MA_C, q = AR_C, measure.type = "CP", measure = list(intensity = "lambda", df = setJumps(jumps = "Gaussian")), XinExpr = XinExpr, Cogarch.var="y", V.var="v", Latent.var="x", ma.par="MA", ar.par="AR")) 
-  }
-  if (name == "Carma Noise: Compound Poisson") {
-    jumpsGlobalEnv <<- jumps
-    model <- yuima::setCarma(p = AR_C, q = MA_C, ma.par="MA", ar.par="AR", measure = list(intensity = "lambda", df = list("dnorm(z, 0, 1)")), scale.par = "sigma_jump", loc.par = "mu_jump", measure.type = "CP", XinExpr = XinExpr)#df = setJumps(jumps = jumpsGlobalEnv)), measure.type = "CP", XinExpr = XinExpr)
-    rm(jumpsGlobalEnv, envir = globalenv())
-    return(model)
-  }
+  if (name == "Cogarch(p,q)") return(yuima::setCogarch(p = MA_C, q = AR_C, measure.type = "CP", measure = list(intensity = "lambda", df = setJumps(jumps = "Gaussian")), XinExpr = XinExpr, Cogarch.var="y", V.var="v", Latent.var="x", ma.par="MA", ar.par="AR")) 
+  if (name == "Carma(p,q)") return(yuima::setCarma(p = AR_C, q = MA_C, ma.par="MA", ar.par="AR", XinExpr = XinExpr))
 }
 
 printModelLatex <- function(names, process, jumps = NA){
@@ -523,7 +517,14 @@
     return (qmle(...))
 }
 
-addModel <- function(modName, modClass, AR_C, MA_C, jumps, symbName, data, delta, start, startMin, startMax, tries, seed, method="BFGS", fixed = list(), lower, upper, joint=FALSE, aggregation=TRUE, threshold=NULL, session, anchorId, alertId){
+clearNA <- function(List){
+  for (i in names(List))
+    if (is.na(List[[i]]))
+      List[[i]] <- NULL
+    return (List)
+}
+
+addModel <- function(modName, modClass, AR_C, MA_C, jumps, symbName, data, delta, start, startMin, startMax, trials, seed, method="BFGS", fixed = list(), lower, upper, joint=FALSE, aggregation=TRUE, threshold=NULL, session, anchorId, alertId){
   info <- list(
     class = modClass,
     modName = modName,
@@ -535,7 +536,7 @@
     start = start,
     startMin = startMin,
     startMax = startMax,
-    tries = tries,
+    trials = trials,
     seed = seed,
     fixed = fixed,
     lower = lower,
@@ -544,12 +545,6 @@
     aggregation = aggregation,
     threshold = threshold
   )
-  clearNA <- function(List){
-    for (i in names(List))
-      if (is.na(List[[i]]))
-        List[[i]] <- NULL
-    return (List)
-  }
   if(!is.na(seed)) set.seed(seed)
   if(is.na(seed)) set.seed(NULL)
   start <- clearNA(start)
@@ -566,29 +561,27 @@
     mu <- alpha +0.5*sigma^2
     if (is.null(start$sigma)) start$sigma <- sigma
     if (is.null(start$mu)) start$mu <- mu
-    QMLE <- try(qmle(model, start = start, fixed = fixed, method = method, lower = lower, upper = upper))
+    QMLE <- try(qmle(model, start = start, fixed = fixed, method = method, lower = lower, upper = upper, rcpp = TRUE))
     if (class(QMLE)=="try-error"){
       createAlert(session = session, anchorId = anchorId, alertId = alertId, content =  paste("Unable to estimate ", modName," on ", symbName, ". Try to use 'Advanced Settings' and customize estimation.", sep = ""), style = "danger")
       return()
     }
   } 
   else if (modClass=="CARMA") {
-    allParam <- unique(c(parameters at drift, parameters at jump, parameters at measure, parameters at xinit[1]))
+    allParam <- unique(c(parameters at drift, parameters at xinit[1]))
     if (all(allParam %in% c(names(start),names(fixed))))
-      QMLE <- try(qmleGUI(model, start = start, method = method, lower = lower, upper = upper, #REMOVE# fixed = fixed, joint = joint, aggregation = aggregation,
-                       threshold = threshold, grideq = TRUE))
+      QMLE <- try(qmleGUI(model, start = start, method = method, lower = lower, upper = upper))
     else {
       miss <- allParam[!(allParam %in% c(names(start),names(fixed)))]
       m2logL_prec <- NA
       na_prec <- NA
       withProgress(message = 'Step: ', value = 0, {
-        for(iter in 1:tries){
-          incProgress(1/tries, detail = paste(iter,"(/", tries ,")"))
+        for(iter in 1:trials){
+          incProgress(1/trials, detail = paste(iter,"(/", trials ,")"))
           for(j in 1:3){
             for (i in miss)
               start[[i]] <- runif(1, min = max(lower[[i]],startMin[[i]], na.rm = TRUE), max = min(upper[[i]],startMax[[i]],na.rm = TRUE))
-            QMLEtemp <- try(qmleGUI(model, start = start, method = method, lower = lower, upper = upper, #fixed = fixed, joint = joint, aggregation = aggregation,
-                                 threshold = threshold, grideq = TRUE))
+            QMLEtemp <- try(qmleGUI(model, start = start, method = method, lower = lower, upper = upper))
             if (class(QMLEtemp)!="try-error") if (all(!is.na(summary(QMLEtemp)@coef[,"Estimate"])))
               break
           }
@@ -598,8 +591,7 @@
               coefTable <- summary(QMLEtemp)@coef
               for (param in rownames(coefTable))
                 start[[param]] <- as.numeric(coefTable[param,"Estimate"])
-              QMLEtemp <- try(qmleGUI(model, start = start, method = method, lower = lower, upper = upper, #fixed = fixed, joint = joint, aggregation = aggregation,
-                               threshold = threshold, grideq = TRUE))
+              QMLEtemp <- try(qmleGUI(model, start = start, method = method, lower = lower, upper = upper))
               if (class(QMLEtemp)=="try-error") break
               else if(summary(QMLEtemp)@m2logL>=m2logL*abs(sign(m2logL)-0.001)) break
             }
@@ -624,7 +616,7 @@
               }
             }
           }
-          if (iter==tries & class(QMLEtemp)=="try-error" & !exists("QMLE")){
+          if (iter==trials & class(QMLEtemp)=="try-error" & !exists("QMLE")){
             createAlert(session = session, anchorId = anchorId, alertId = alertId, content =  paste("Unable to estimate ", modName," on ", symbName, ". Try to use 'Advanced Settings' and customize estimation.", sep = ""), style = "danger")
             return()
           }
@@ -636,19 +628,19 @@
     allParam <- unique(c(parameters at drift, parameters at xinit))
     if (all(allParam %in% c(names(start),names(fixed))))
       QMLE <- try(qmle(model, start = start, fixed = fixed, method = method, lower = lower, upper = upper, #REMOVE# joint = joint, aggregation = aggregation,
-                       threshold = threshold, grideq = TRUE))
+                       threshold = threshold, grideq = TRUE, rcpp = TRUE))
     else {
       miss <- allParam[!(allParam %in% c(names(start),names(fixed)))]
       m2logL_prec <- NA
       na_prec <- NA
       withProgress(message = 'Step: ', value = 0, {
-        for(iter in 1:tries){
-          incProgress(1/tries, detail = paste(iter,"(/", tries ,")"))
+        for(iter in 1:trials){
+          incProgress(1/trials, detail = paste(iter,"(/", trials ,")"))
           for(j in 1:3){
             for (i in miss)
               start[[i]] <- runif(1, min = max(lower[[i]],startMin[[i]], na.rm = TRUE), max = min(upper[[i]],startMax[[i]],na.rm = TRUE))
             QMLEtemp <- try(qmle(model, start = start, fixed = fixed, method = method, lower = lower, upper = upper, #joint = joint, aggregation = aggregation,
-                                 threshold = threshold, grideq = TRUE))
+                                 threshold = threshold, grideq = TRUE, rcpp = TRUE))
             if (class(QMLEtemp)!="try-error") if (all(!is.na(summary(QMLEtemp)@coef[,"Estimate"])))
               break
           }
@@ -659,7 +651,7 @@
               for (param in rownames(coefTable))
                 start[[param]] <- as.numeric(coefTable[param,"Estimate"])
               QMLEtemp <- try(qmle(model, start = start, fixed = fixed, method = method, lower = lower, upper = upper, #joint = joint, aggregation = aggregation,
-                                   threshold = threshold, grideq = TRUE))
+                                   threshold = threshold, grideq = TRUE, rcpp = TRUE))
               if (class(QMLEtemp)=="try-error") break
               else if(summary(QMLEtemp)@objFunVal>=m2logL*abs(sign(m2logL)-0.001)) break
             }
@@ -684,7 +676,7 @@
               }
             }
           }
-          if (iter==tries & class(QMLEtemp)=="try-error" & !exists("QMLE")){
+          if (iter==trials & class(QMLEtemp)=="try-error" & !exists("QMLE")){
             createAlert(session = session, anchorId = anchorId, alertId = alertId, content =  paste("Unable to estimate ", modName," on ", symbName, ". Try to use 'Advanced Settings' and customize estimation.", sep = ""), style = "danger")
             return()
           }
@@ -703,8 +695,8 @@
       m2logL_prec <- NA
       na_prec <- NA
       withProgress(message = 'Step: ', value = 0, {
-        for(iter in 1:tries){
-          incProgress(1/tries, detail = paste(iter,"(/", tries ,")"))
+        for(iter in 1:trials){
+          incProgress(1/trials, detail = paste(iter,"(/", trials ,")"))
           for(j in 1:3){
             for (i in miss)
               start[[i]] <- runif(1, min = max(lower[[i]],startMin[[i]], na.rm = TRUE), max = min(upper[[i]],startMax[[i]],na.rm = TRUE))
@@ -745,7 +737,7 @@
               }
             }
           }
-          if (iter==tries & class(QMLEtemp)=="try-error" & !exists("QMLE")){
+          if (iter==trials & class(QMLEtemp)=="try-error" & !exists("QMLE")){
             createAlert(session = session, anchorId = anchorId, alertId = alertId, content =  paste("Unable to estimate ", modName," on ", symbName, ". Try to use 'Advanced Settings' and customize estimation.", sep = ""), style = "danger")
             return()
           }
@@ -756,19 +748,19 @@
   else {
     if (all(parameters at all %in% c(names(start),names(fixed))))
       QMLE <- try(qmle(model, start = start, fixed = fixed, method = method, lower = lower, upper = upper, #REMOVE# joint = joint, aggregation = aggregation,
-                       threshold = threshold))
+                       threshold = threshold, rcpp = TRUE))
     else {
       miss <- parameters at all[!(parameters at all %in% c(names(start),names(fixed)))]
       m2logL_prec <- NA
       na_prec <- NA
       withProgress(message = 'Step: ', value = 0, {
-        for(iter in 1:tries){
-          incProgress(1/tries, detail = paste(iter,"(/", tries ,")"))
+        for(iter in 1:trials){
+          incProgress(1/trials, detail = paste(iter,"(/", trials ,")"))
           for(j in 1:3){
             for (i in miss)
               start[[i]] <- runif(1, min = max(lower[[i]],startMin[[i]], na.rm = TRUE), max = min(upper[[i]],startMax[[i]],na.rm = TRUE))
             QMLEtemp <- try(qmle(model, start = start, fixed = fixed, method = method, lower = lower, upper = upper, #joint = joint, aggregation = aggregation,
-                                 threshold = threshold))
+                                 threshold = threshold, rcpp = TRUE))
             if (class(QMLEtemp)!="try-error") if (all(!is.na(summary(QMLEtemp)@coef[,"Estimate"])))
               break
           }
@@ -779,7 +771,7 @@
               for (param in names(start))
                 start[[param]] <- as.numeric(coefTable[param,"Estimate"])
               QMLEtemp <- try(qmle(model, start = start, fixed = fixed, method = method, lower = lower, upper = upper, #joint = joint, aggregation = aggregation,
-                                   threshold = threshold))
+                                   threshold = threshold, rcpp = TRUE))
               if (class(QMLEtemp)=="try-error") break
               else if (summary(QMLEtemp)@m2logL>=m2logL*abs(sign(m2logL)-0.001)) break
             }
@@ -804,7 +796,7 @@
               }
             }
           }
-          if (iter==tries & class(QMLEtemp)=="try-error" & !exists("QMLE")){
+          if (iter==trials & class(QMLEtemp)=="try-error" & !exists("QMLE")){
             createAlert(session = session, anchorId = anchorId, alertId = alertId, content =  paste("Unable to estimate ", modName," on ", symbName, ". Try to use 'Advanced Settings' and customize estimation.", sep = ""), style = "danger")
             return()
           }
@@ -813,10 +805,8 @@
     }
   }
   
+  if(!exists("QMLE")) return()
   
-  if(!exists("QMLE"))
-    return()
-  
   yuimaGUIdata$model[[symbName]][[ifelse(is.null(length(yuimaGUIdata$model[[symbName]])),1,length(yuimaGUIdata$model[[symbName]])+1)]] <<- list(
    model = model,
    qmle = QMLE,
@@ -828,8 +818,137 @@
 
 
 
+addCPoint <- function(modelName, symb, trials, frac = 0.15, delta = 0.01, session, anchorId, alertId = NULL){
+  series <- getData(symb)
+  mod <- setModelByName(name = modelName)
+  bounds <- defaultBounds(name = modelName)
+  startBounds <- defaultBounds(name = modelName, lower = -100, upper = 100)
+  yuima <- setYuima(data = setData(series, delta = delta), model = mod)
+  start <- list()
+  startMin <- startBounds$lower
+  startMax <- startBounds$upper
+  lower <- clearNA(bounds$lower)
+  upper <- clearNA(bounds$upper)
+  miss <- mod at parameter@all
+  
+  m2logL_prec <- NA
+  na_prec <- NA
+  for(iter in 1:trials){
+    for(j in 1:3){
+      for (i in miss)
+        start[[i]] <- runif(1, min = max(lower[[i]],startMin[[i]], na.rm = TRUE), max = min(upper[[i]],startMax[[i]],na.rm = TRUE))
+      QMLEtempL <- try(qmleL(yuima = yuima, t = frac*length(series)*delta, start = start, method="L-BFGS-B", lower = lower, upper = upper, rcpp = TRUE))
+      if (class(QMLEtempL)!="try-error") if (all(!is.na(summary(QMLEtempL)@coef[,"Estimate"])))
+        break
+    }
+    if (class(QMLEtempL)!="try-error") if (all(!is.na(summary(QMLEtempL)@coef[,"Estimate"]))){
+      repeat{
+        m2logL <- summary(QMLEtempL)@m2logL
+        coefTable <- summary(QMLEtempL)@coef
+        for (param in names(start))
+          start[[param]] <- as.numeric(coefTable[param,"Estimate"])
+        QMLEtempL <- try(qmleL(yuima = yuima, t = frac*length(series)*delta, start = start, method="L-BFGS-B", lower = lower, upper = upper, rcpp = TRUE))
+        if (class(QMLEtempL)=="try-error") break
+        else if (summary(QMLEtempL)@m2logL>=m2logL*abs(sign(m2logL)-0.001)) break
+      }
+      if(is.na(m2logL_prec) & class(QMLEtempL)!="try-error"){
+        QMLEL <- QMLEtempL
+        m2logL_prec <- summary(QMLEL)@m2logL
+        na_prec <- sum(is.na(coefTable))
+      }
+      else if (class(QMLEtempL)!="try-error"){
+        if (sum(is.na(coefTable)) < na_prec){
+          QMLEL <- QMLEtempL
+          m2logL_prec <- summary(QMLEL)@m2logL
+          na_prec <- sum(is.na(coefTable))
+        }
+        else {
+          test <- summary(QMLEtempL)@m2logL
+          if(test < m2logL_prec & sum(is.na(coefTable))==na_prec){
+            QMLEL <- QMLEtempL
+            m2logL_prec <- test
+            na_prec <- sum(is.na(coefTable))
+          }
+        }
+      }
+    }
+  }
+  if (iter==trials & class(QMLEtempL)=="try-error" & !exists("QMLEL")){
+    createAlert(session = session, anchorId = anchorId, alertId = alertId, content =  paste("Unable to estimate change points of ", symb, ". Try to increase the number of Trials", sep = ""), style = "error")
+    return()
+  }
+  if(!exists("QMLEL")) return()
+  tmpL <- QMLEL
 
+  m2logL_prec <- NA
+  na_prec <- NA
+  for(iter in 1:trials){
+    for(j in 1:3){
+      for (i in miss)
+        start[[i]] <- runif(1, min = max(lower[[i]],startMin[[i]], na.rm = TRUE), max = min(upper[[i]],startMax[[i]],na.rm = TRUE))
+      QMLEtempR <- try(qmleR(yuima = yuima, t = (1-frac)*length(series)*delta, start = start, method="L-BFGS-B", lower = lower, upper = upper, rcpp = TRUE))
+      if (class(QMLEtempR)!="try-error") if (all(!is.na(summary(QMLEtempR)@coef[,"Estimate"])))
+        break
+    }
+    if (class(QMLEtempR)!="try-error") if (all(!is.na(summary(QMLEtempR)@coef[,"Estimate"]))){
+      repeat{
+        m2logL <- summary(QMLEtempR)@m2logL
+        coefTable <- summary(QMLEtempR)@coef
+        for (param in names(start))
+          start[[param]] <- as.numeric(coefTable[param,"Estimate"])
+        QMLEtempR <- try(qmleR(yuima = yuima, t = (1-frac)*length(series)*delta, start = start, method="L-BFGS-B", lower = lower, upper = upper, rcpp = TRUE))
+        if (class(QMLEtempR)=="try-error") break
+        else if (summary(QMLEtempR)@m2logL>=m2logL*abs(sign(m2logL)-0.001)) break
+      }
+      if(is.na(m2logL_prec) & class(QMLEtempR)!="try-error"){
+        QMLER <- QMLEtempR
+        m2logL_prec <- summary(QMLER)@m2logL
+        na_prec <- sum(is.na(coefTable))
+      }
+      else if (class(QMLEtempR)!="try-error"){
+        if (sum(is.na(coefTable)) < na_prec){
+          QMLER <- QMLEtempR
+          m2logL_prec <- summary(QMLER)@m2logL
+          na_prec <- sum(is.na(coefTable))
+        }
+        else {
+          test <- summary(QMLEtempR)@m2logL
+          if(test < m2logL_prec & sum(is.na(coefTable))==na_prec){
+            QMLER <- QMLEtempR
+            m2logL_prec <- test
+            na_prec <- sum(is.na(coefTable))
+          }
+        }
+      }
+    }
+  }
+  if (iter==trials & class(QMLEtempR)=="try-error" & !exists("QMLER")){
+    createAlert(session = session, anchorId = anchorId, alertId = alertId, content =  paste("Unable to estimate change points of ", symb, ". Try to increase the number of Trials", sep = ""), style = "error")
+    return()
+  }
+  if(!exists("QMLER")) return()
+  tmpR <- QMLER
+  
+  cp_prec <- CPoint(yuima = yuima, param1=coef(tmpL), param2=coef(tmpR))
+  repeat{
+    tmpL <- qmleL(yuima, start=as.list(coef(tmpL)), t = cp_prec$tau, lower=lower, upper = upper, method="L-BFGS-B", rcpp = TRUE)
+    tmpR <- qmleR(yuima, start=as.list(coef(tmpR)), t = cp_prec$tau, lower=lower, upper = upper, method="L-BFGS-B", rcpp = TRUE)
+    cp <- CPoint(yuima = yuima, param1=coef(tmpL), param2=coef(tmpR))
+    if (abs(cp$tau - cp_prec$tau)<delta) break
+    else cp_prec <- cp
+  }
+  
+  yuimaGUIdata$cpYuima[[symb]] <<- list(tau = index(series)[as.integer(cp$tau/delta)], model = modelName, trials = trials)
+  
+}
 
+
+
+
+
+
+
+
 getModelNames <- function(){
   return(isolate({yuimaGUItable$model}))
 }
@@ -861,13 +980,16 @@
   }
   is.valid <- TRUE
   model <- modelYuima at model
-  if (info$class=="COGARCH") incr.L <- cogarchNoise(yuima = modelYuima, param = true.parameter)$incr.L
+  if (info$class=="COGARCH") increments <- cogarchNoise(yuima = modelYuima, param = true.parameter)$incr.L
+  if (info$class=="CARMA") increments <- CarmaNoise(yuima = modelYuima, param = true.parameter)
   withProgress(message = 'Simulating: ', value = 0, {
     for (i in 1:nsim){
       incProgress(1/nsim, detail = paste("Simulating:",i,"(/",nsim,")"))
       if (info$class=="COGARCH")
-        simulation <- try(yuima::simulate(object = model, increment.L = sample(x = incr.L, size = sampling at n, replace = TRUE), xinit = xinit, true.parameter = true.parameter, nsim = nsim, sampling = sampling, space.discretized = space.discretized, method = method))
-      else 
+        simulation <- try(yuima::simulate(object = model, increment.L = sample(x = increments, size = sampling at n, replace = TRUE), xinit = xinit, true.parameter = true.parameter, nsim = nsim, sampling = sampling, space.discretized = space.discretized, method = method))
+      else if (info$class=="CARMA")
+        simulation <- try(yuima::simulate(object = model, increment.W = t(sample(x = increments, size = sampling at n, replace = TRUE)), xinit = xinit, true.parameter = true.parameter, nsim = nsim, sampling = sampling, space.discretized = space.discretized, method = method))
+      else
         simulation <- try(yuima::simulate(object = model, xinit = xinit, true.parameter = true.parameter, nsim = nsim, sampling = sampling, space.discretized = space.discretized, method = method))
       if (class(simulation)=="try-error"){
         is.valid <- FALSE
@@ -887,10 +1009,11 @@
   }
 
   if(saveTraj==TRUE){
+    times <- index(trajectory)
     if(class(info$simulate.from)=="Date")
-      index(trajectory) <- as.POSIXct(24*60*60*index(trajectory)/simulation at sampling@delta*as.numeric(info$simulate.to-info$simulate.from)/(simulation at sampling@n), origin = info$simulate.from)
+      index(trajectory) <- as.POSIXct(24*60*60*(times-times[1])/simulation at sampling@delta*as.numeric(info$simulate.to-info$simulate.from)/(simulation at sampling@n), origin = info$simulate.from)
     if(class(info$simulate.from)=="numeric")
-      index(trajectory) <- as.numeric(index(trajectory)/simulation at sampling@delta*as.numeric(info$simulate.to-info$simulate.from)/(simulation at sampling@n))
+      index(trajectory) <- as.numeric(times/simulation at sampling@delta*as.numeric(info$simulate.to-info$simulate.from)/(simulation at sampling@n))
     if(!is.null(colnames(trajectory)))
       colnames(trajectory) <- seq(1:length(colnames(trajectory)))
   }
@@ -939,6 +1062,9 @@
   }
 }
 
+
+
+
 addHedging <- function(model, symbName, info = list(), xinit, true.parameter, nsim, sampling, session, anchorId){
   closeAlert(session, "addHedging_alert")
   hist <- vector()
@@ -1048,15 +1174,11 @@
 
 
 
-CPanalysis <- function(x, method = c("lSQ", "KSdiff", "KSperc"), pvalue = 0.01){
+CPanalysis <- function(x, method = c("KSdiff", "KSperc"), pvalue = 0.01){
   if (pvalue > 0.1){
     pvalue <- 0.1
     warning("pvalue re-defined: 0.1")
   }
-  if(method=="lSQ"){
-    tau <- cpoint(x)$tau0
-    return(list(tau=tau, pvalue=NA))
-  }
   if(method=="KSdiff" | method=="KSperc"){
     x_incr <- switch (method,
                       "KSdiff" = na.omit(diff(x)),
@@ -1087,7 +1209,7 @@
       tau <- NA
       p.value <- NA
     }
-    return (list(tau=tau,pvalue=p.value))
+    return (list(tau=tau,pvalue=p.value, method=method))
   }  
 }
 

Modified: pkg/yuimaGUI/inst/yuimaGUI/server.R
===================================================================
--- pkg/yuimaGUI/inst/yuimaGUI/server.R	2016-10-15 15:14:36 UTC (rev 482)
+++ pkg/yuimaGUI/inst/yuimaGUI/server.R	2016-10-16 19:05:53 UTC (rev 483)
@@ -255,7 +255,7 @@
   })
   
   output$jumps <- renderUI({
-    if (input$modelClass!="Diffusion process" & input$modelClass!="COGARCH")
+    if (input$modelClass!="Diffusion process" & input$modelClass!="COGARCH" & input$modelClass!="CARMA")
       return(selectInput("jumps",label = "Jumps", choices = defaultJumps))
   })
   
@@ -460,38 +460,72 @@
   ###Choose Range input set to "Select range from charts" if charts have been brushed
   output$chooseRange <- renderUI({
     sel <- "full"
-    if (!is.null(range_selectRange$x))
-      sel <- "selected"
-    selectInput("chooseRange", label = "Range", choices = c("Full Range" = "full", "Select range from charts" = "selected"), selected = sel)
+    if (!is.null(range_selectRange$x)) sel <- "selected"
+    selectInput("chooseRange", label = "Range", choices = c("Full Range" = "full", "Select Range from Charts" = "selected", "Specify Range" = "specify"), selected = sel)
   })
+  
+  output$chooseRange_specify <- renderUI({
+    if(!is.null(input$plotsRangeSeries)) {
+      data <- getData(input$plotsRangeSeries)
+    if(class(index(data))=="numeric") 
+      return(div(
+        column(6,numericInput("chooseRange_specify_t0", label = "From", min = start(data), max = end(data), value = start(data))),
+        column(6,numericInput("chooseRange_specify_t1", label = "To", min = start(data), max = end(data), value = end(data)))
+      ))
+    if(class(index(data))=="Date")
+      return(dateRangeInput("chooseRange_specify_date", start = start(data), end = end(data), label = "Specify Range"))
+    }
+  })
 
+  
+  observe({
+    shinyjs::toggle(id = "chooseRange_specify", condition = (input$chooseRange)=="specify")
+  })
+
   ###Function to update data range to use to estimate models
-  updateRange_seriesToEstimate <- function(symb, range = c("full","selected"), type = c("Date", "numeric")){
+  updateRange_seriesToEstimate <- function(symb, range = c("full","selected","specify"), type = c("Date", "numeric")){
     for (i in symb){
-        data <- getData(i)
-        if (range == "full"){
-          levels(seriesToEstimate$table[,"From"]) <- c(levels(seriesToEstimate$table[,"From"]), as.character(start(data)))
-          levels(seriesToEstimate$table[,"To"]) <- c(levels(seriesToEstimate$table[,"To"]), as.character(end(data)))
-          seriesToEstimate$table[i,"From"] <<- as.character(start(data))
-          seriesToEstimate$table[i,"To"] <<- as.character(end(data))
+      data <- getData(i)
+      if (range == "full"){
+        levels(seriesToEstimate$table[,"From"]) <- c(levels(seriesToEstimate$table[,"From"]), as.character(start(data)))
+        levels(seriesToEstimate$table[,"To"]) <- c(levels(seriesToEstimate$table[,"To"]), as.character(end(data)))
+        seriesToEstimate$table[i,"From"] <<- as.character(start(data))
+        seriesToEstimate$table[i,"To"] <<- as.character(end(data))
+      }
+      if (range == "selected"){
+        if(!is.null(range_selectRange$x) & class(index(data))==type){
+          start <- range_selectRange$x[1]
+          end <- range_selectRange$x[2]
+          if(class(index(data))=="numeric"){
+            start <- as.numeric(start)
+            end <- as.numeric(end)
+          }
+          start <- max(start(data),start)
+          end <- min(end(data), end)
+          levels(seriesToEstimate$table[,"From"]) <- c(levels(seriesToEstimate$table[,"From"]), as.character(start))
+          levels(seriesToEstimate$table[,"To"]) <- c(levels(seriesToEstimate$table[,"To"]), as.character(end))
+          seriesToEstimate$table[i,"From"] <<- as.character(start)
+          seriesToEstimate$table[i,"To"] <<- as.character(end)
         }
-        if (range == "selected"){
-          if(!is.null(range_selectRange$x) & class(index(data))==type){
-            start <- range_selectRange$x[1]
-            end <- range_selectRange$x[2]
-            if(class(index(data))=="numeric"){
-              start <- as.integer(start)
-              end <- as.integer(end)
-            }
-            start <- max(start(data),start)
-            end <- min(end(data), end)
-            levels(seriesToEstimate$table[,"From"]) <- c(levels(seriesToEstimate$table[,"From"]), as.character(start))
-            levels(seriesToEstimate$table[,"To"]) <- c(levels(seriesToEstimate$table[,"To"]), as.character(end))
-            seriesToEstimate$table[i,"From"] <<- as.character(start)
-            seriesToEstimate$table[i,"To"] <<- as.character(end)
+      }
+      if (range == "specify"){
+        if(class(index(data))==type){
+          if(class(index(data))=="Date"){
+            start <- input$chooseRange_specify_date[1]
+            end <- input$chooseRange_specify_date[2]
           }
+          if(class(index(data))=="numeric"){
+            start <- input$chooseRange_specify_t0
+            end <- input$chooseRange_specify_t1
+          }
+          start <- max(start(data),start)
+          end <- min(end(data), end)
+          levels(seriesToEstimate$table[,"From"]) <- c(levels(seriesToEstimate$table[,"From"]), as.character(start))
+          levels(seriesToEstimate$table[,"To"]) <- c(levels(seriesToEstimate$table[,"To"]), as.character(end))
+          seriesToEstimate$table[i,"From"] <<- as.character(start)
+          seriesToEstimate$table[i,"To"] <<- as.character(end)
         }
-
+      }
     }
   }
 
@@ -517,7 +551,7 @@
       if (is.null(deltaSettings[[symb]]))
         deltaSettings[[symb]] <<- 0.01
       for (modName in input$model){
-        if (class(try(setModelByName(modName, jumps = switch(input$modelClass, "Diffusion process" = NA, "Compound Poisson" = input$jumps, "COGARCH"=NA, "CARMA"=input$jumps), AR_C = ifelse(input$modelClass %in% c("CARMA","COGARCH"), input$AR_C, NA), MA_C = ifelse(input$modelClass %in% c("CARMA","COGARCH"), input$MA_C, NA))))!="try-error"){
+        if (class(try(setModelByName(modName, jumps = switch(input$modelClass, "Diffusion process" = NA, "Compound Poisson" = input$jumps, "COGARCH"=NA, "CARMA"=NA), AR_C = ifelse(input$modelClass %in% c("CARMA","COGARCH"), input$AR_C, NA), MA_C = ifelse(input$modelClass %in% c("CARMA","COGARCH"), input$MA_C, NA))))!="try-error"){
           if (is.null(estimateSettings[[modName]]))
             estimateSettings[[modName]] <<- list()
           if (is.null(estimateSettings[[modName]][[symb]]))
@@ -528,34 +562,34 @@
             estimateSettings[[modName]][[symb]][["start"]] <<- list()
           if (is.null(estimateSettings[[modName]][[symb]][["startMin"]]) | isolate({input$modelClass!="Diffusion process"}))
[TRUNCATED]

To get the complete diff run:
    svnlook diff /svnroot/yuima -r 483


More information about the Yuima-commits mailing list