[Yuima-commits] r751 - pkg/yuima/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Mar 15 06:02:13 CET 2021
Author: yumauehara
Date: 2021-03-15 06:02:13 +0100 (Mon, 15 Mar 2021)
New Revision: 751
Modified:
pkg/yuima/R/qmle.R
Log:
fixed a bug
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2021-03-13 02:39:54 UTC (rev 750)
+++ pkg/yuima/R/qmle.R 2021-03-15 05:02:13 UTC (rev 751)
@@ -16,16 +16,16 @@
DRIFT <- yuima at model@drift
# n <- length(yuima)[1]
n <- dim(env$X)[1]
-
+
drift <- matrix(0,n,d.size)
tmp.env <- new.env()
assign(yuima at model@time.variable, env$time, envir=tmp.env)
-
-
+
+
for(i in 1:length(theta)){
assign(names(theta)[i],theta[[i]], envir=tmp.env)
}
-
+
for(d in 1:d.size){
assign(modelstate[d], env$X[,d], envir=tmp.env)
}
@@ -32,7 +32,7 @@
for(d in 1:d.size){
drift[,d] <- eval(DRIFT[d], envir=tmp.env)
}
-
+
return(drift)
}
@@ -50,11 +50,11 @@
for(i in 1:length(theta)){
assign(names(theta)[i],theta[[i]],envir=tmp.env)
}
-
+
for(d in 1:d.size){
assign(modelstate[d], env$X[,d], envir=tmp.env)
}
-
+
for(r in 1:r.size){
for(d in 1:d.size){
diff[d, r, ] <- eval(DIFFUSION[[d]][r], envir=tmp.env)
@@ -72,7 +72,7 @@
d.size <- yuima at model@equation.number
modelstate <- yuima at model@state.variable
n <- dim(env$X)[1]
-
+
tmp.env <- new.env()
assign(yuima at model@time.variable, env$time, envir =tmp.env)
JUMP <- yuima at model@jump.coeff
@@ -80,7 +80,7 @@
for(i in 1:length(theta)){
assign(names(theta)[i],theta[[i]],envir=tmp.env)
}
-
+
for(d in 1:d.size){
assign(modelstate[d], env$X[,d],envir=tmp.env)
}
@@ -140,16 +140,16 @@
method<-"L-BFGS-B"
}
call <- match.call()
-
+
if( missing(yuima))
yuima.stop("yuima object is missing.")
if(is.COGARCH(yuima)){
if(missing(lower))
lower <- list()
-
+
if(missing(upper))
upper <- list()
-
+
res <- NULL
if("grideq" %in% names(as.list(call)[-(1:2)])){
res <- PseudoLogLik.COGARCH(yuima, start, method=method, fixed = list(),
@@ -158,17 +158,17 @@
res <- PseudoLogLik.COGARCH(yuima, start, method=method, fixed = list(),
lower, upper, Est.Incr, call, grideq = FALSE, aggregation = aggregation,...)
}
-
+
return(res)
}
-
+
if(is.PPR(yuima)){
if(missing(lower))
lower <- list()
-
+
if(missing(upper))
upper <- list()
-
+
# res <- NULL
# if("grideq" %in% names(as.list(call)[-(1:2)])){
res <- quasiLogLik.PPR(yuimaPPR = yuima, parLambda = start, method=method, fixed = list(),
@@ -177,43 +177,43 @@
# res <- PseudoLogLik.COGARCH(yuima, start, method=method, fixed = list(),
# lower, upper, Est.Incr, call, grideq = FALSE, aggregation = aggregation,...)
# }
-
+
return(res)
}
-
+
orig.fixed <- fixed
orig.fixed.par <- names(orig.fixed)
if(is.Poisson(yuima))
threshold <- 0
## param handling
-
+
## FIXME: maybe we should choose initial values at random within lower/upper
## at present, qmle stops
if( missing(start) )
yuima.stop("Starting values for the parameters are missing.")
-
+
#14/12/2013 We modify the QMLE function when the model is a Carma(p,q).
# In this case we use a two step procedure:
# First) The Coefficient are obtained by QMLE computed using the Kalman Filter.
# Second) Using the result in Brockwell, Davis and Yang (2007) we retrieve
# the underlying Levy. The estimated increments are used to find the L?vy parameters.
-
+
# if(is(yuima at model, "yuima.carma")){
# yuima.warm("two step procedure for carma(p,q)")
# return(null)
# }
#
-
+
yuima.nobs <- as.integer(max(unlist(lapply(get.zoo.data(yuima),length))-1,na.rm=TRUE))
-
+
diff.par <- yuima at model@parameter at diffusion
-
+
# 24/12
if(is.CARMA(yuima) && length(diff.par)==0
&& length(yuima at model@parameter at jump)!=0){
diff.par<-yuima at model@parameter at jump
}
-
+
if(is.CARMA(yuima) && length(yuima at model@parameter at jump)!=0){
CPlist <- c("dgamma", "dexp")
codelist <- c("rIG", "rgamma")
@@ -229,17 +229,17 @@
}
# return(NULL)
}
-
+
}
-
+
if(yuima at model@measure.type=="code"){
if(class(yuima at model@measure$df)=="yuima.law"){
measurefunc <- "yuima.law"
}
else{
-
- tmp <- regexpr("\\(", yuima at model@measure$df$exp)[1]
- measurefunc <- substring(yuima at model@measure$df$exp, 1, tmp-1)
+
+ tmp <- regexpr("\\(", yuima at model@measure$df$exp)[1]
+ measurefunc <- substring(yuima at model@measure$df$exp, 1, tmp-1)
}
if(!is.na(match(measurefunc,codelist))){
yuima.warn("carma(p,q): the qmle for a carma(p,q) driven by a non-Negative Levy will be implemented as soon as possible")
@@ -250,29 +250,29 @@
#return(NULL)
}
}
-
-
+
+
# yuima.warn("carma(p,q): the qmle for a carma(p,q) driven by a Jump process will be implemented as soon as possible ")
# return(NULL)
}
-
+
# 24/12
if(is.CARMA(yuima) && length(yuima at model@info at lin.par)>0){
yuima.warn("carma(p,q): the case of lin.par will be implemented as soon as")
return(NULL)
}
-
-
+
+
drift.par <- yuima at model@parameter at drift
#01/01 we introduce the new variable in order
# to take into account the parameters in the starting conditions
-
+
if(is.CARMA(yuima)){
#if(length(yuima at model@info at scale.par)!=0){
xinit.par <- yuima at model@parameter at xinit
#}
}
-
+
# SMI-2/9/14: measure.par is used for Compound Poisson
# and CARMA, jump.par only by CARMA
jump.par <- NULL
@@ -288,7 +288,7 @@
}
# jump.par is used for CARMA
common.par <- yuima at model@parameter at common
-
+
JointOptim <- joint
if(is.CARMA(yuima) && length(yuima at model@parameter at jump)!=0){
if(any((match(jump.par, drift.par)))){
@@ -296,8 +296,8 @@
yuima.warn("Drift and diffusion parameters must be different. Doing
joint estimation, asymptotic theory may not hold true.")
}
- }
-
+ }
+
if(length(common.par)>0){
JointOptim <- TRUE
yuima.warn("Drift and diffusion parameters must be different. Doing
@@ -309,32 +309,32 @@
# #return(NULL)
# }
}
-
+
# if(!is(yuima at model, "yuima.carma")){
- # if(length(jump.par)+length(measure.par)>0)
+ # if(length(jump.par)+yuima at model@measure.type == "CP")
# yuima.stop("Cannot estimate the jump models, yet")
# }
-
-
+
+
if(!is.list(start))
yuima.stop("Argument 'start' must be of list type.")
-
+
fullcoef <- NULL
-
+
if(length(diff.par)>0)
fullcoef <- diff.par
-
+
if(length(drift.par)>0)
fullcoef <- c(fullcoef, drift.par)
-
+
if(is.CARMA(yuima) &&
(length(yuima at model@info at loc.par)!=0)){
# 01/01 We modify the code for considering
# the loc.par in yuima.carma model
fullcoef<-c(fullcoef, yuima at model@info at loc.par)
-
+
}
-
+
if(is.CARMA(yuima) && (NoNeg.Noise==TRUE)){
if((yuima at model@info at q+1)==yuima at model@info at p){
mean.noise<-"mean.noise"
@@ -341,11 +341,11 @@
fullcoef<-c(fullcoef, mean.noise)
}
}
-
- # if(is.CARMA(yuima) && (length(measure.par)>0)){
+
+ # if(is.CARMA(yuima) && (yuima at model@measure.type == "CP")){
fullcoef<-c(fullcoef, measure.par)
#}
-
+
if(is.CARMA(yuima)){
if(length(yuima at model@parameter at xinit)>1){
#fullcoef<-unique(c(fullcoef,yuima at model@parameter at xinit))
@@ -356,14 +356,14 @@
}
}
}
-
-
+
+
npar <- length(fullcoef)
-
-
+
+
fixed.par <- names(fixed) # We use Fixed.par when we consider a Carma with scale parameter
fixed.carma=NULL
- if(is.CARMA(yuima) && (length(measure.par)>0)){
+ if(is.CARMA(yuima) && (length(measure.par) > 0)){
if(!missing(fixed)){
if(names(fixed) %in% measure.par){
idx.fixed.carma<-match(names(fixed),measure.par)
@@ -396,10 +396,10 @@
}
}
}
-
-
-
-
+
+
+
+
for( j in c(1:length(measure.par))){
if(is.na(match(measure.par[j],names(fixed)))){
fixed.par <- c(fixed.par,measure.par[j])
@@ -406,21 +406,21 @@
fixed[measure.par[j]]<-start[measure.par[j]]
}
}
-
+
}
if (any(!(fixed.par %in% fullcoef)))
yuima.stop("Some named arguments in 'fixed' are not arguments to the supplied yuima model")
-
+
nm <- names(start)
-
+
oo <- match(nm, fullcoef)
-
+
if(any(is.na(oo)))
yuima.stop("some named arguments in 'start' are not arguments to the supplied yuima model")
start <- start[order(oo)]
nm <- names(start)
-
-
+
+
idx.diff <- match(diff.par, nm)
idx.drift <- match(drift.par, nm)
# SMI-2/9/14: idx.measure for CP
@@ -432,13 +432,13 @@
#}
}
#if(is.null(fixed.carma)){
- idx.fixed <- match(fixed.par, nm)
-# }else{
+ idx.fixed <- match(fixed.par, nm)
+ # }else{
# dummynm <- nm[!(nm %in% fixed.par)]
# idx.fixed <- match(fixed.par, dummynm)
# }
orig.idx.fixed <- idx.fixed
-
+
tmplower <- as.list( rep( -Inf, length(nm)))
names(tmplower) <- nm
if(!missing(lower)){
@@ -448,7 +448,7 @@
tmplower[ idx ] <- lower
}
lower <- tmplower
-
+
tmpupper <- as.list( rep( Inf, length(nm)))
names(tmpupper) <- nm
if(!missing(upper)){
@@ -458,11 +458,11 @@
tmpupper[ idx ] <- upper
}
upper <- tmpupper
-
-
-
-
-
+
+
+
+
+
d.size <- yuima at model@equation.number
if (is.CARMA(yuima)){
# 24/12
@@ -469,16 +469,16 @@
d.size <-1
}
n <- length(yuima)[1]
-
+
env <- new.env()
-
+
assign("X", as.matrix(onezoo(yuima)), envir=env)
assign("deltaX", matrix(0, n-1, d.size), envir=env)
# SMI-2/9/14: for CP
assign("Cn.r", numeric(n-1), envir=env)
- if(length(measure.par)==0)
+ if(length(yuima at model@measure.type) == 0)
threshold <- 0 # there are no jumps, we take all observations
-
+
if (is.CARMA(yuima)){
#24/12 If we consider a carma model,
# the observations are only the first column of env$X
@@ -494,8 +494,8 @@
assign("p", yuima at model@info at p, envir=env)
assign("q", yuima at model@info at q, envir=env)
assign("V_inf0", matrix(diag(rep(1,env$p)),env$p,env$p), envir=env)
-
-
+
+
# env$X<-as.matrix(env$X[,1])
# env$deltaX<-as.matrix(env$deltaX[,1])
# assign("time.obs",length(env$X), envir=env)
@@ -503,23 +503,23 @@
# assign("V_inf0", matrix(diag(rep(1,p)),p,p), envir=env)
}
assign("time", as.numeric(index(yuima at data@zoo.data[[1]])), envir=env)
-
+
for(t in 1:(n-1)){
env$deltaX[t,] <- env$X[t+1,] - env$X[t,]
if(!is.CARMA(yuima))
env$Cn.r[t] <- ((sqrt( env$deltaX[t,] %*% env$deltaX[t,])) <= threshold)
}
-
- if(length(measure.par)==0)
+
+ if(length(yuima at model@measure.type) == 0)
env$Cn.r <- rep(1, length(env$Cn.r)) # there are no jumps, we take all observations
-
+
assign("h", deltat(yuima at data@zoo.data[[1]]), envir=env)
-
+
#SMI: 2/9/214 jump
- if(length(measure.par)>0){
-
- # "yuima.law" LM 13/05/2018
+ if(length(yuima at model@measure.type) > 0 && yuima at model@measure.type == "CP"){
+ # "yuima.law" LM 13/05/2018
+
if(class(yuima at model@measure$df)=="yuima.law"){
args <- yuima at model@parameter at measure
}else{
@@ -526,15 +526,17 @@
args <- unlist(strsplit(suppressWarnings(sub("^.+?\\((.+)\\)", "\\1",yuima at model@measure$df$expr,perl=TRUE)), ","))
}
idx.intensity <- numeric(0)
+ if(length(measure.par) > 0){
for(i in 1:length(measure.par)){
if(sum(grepl(measure.par[i],yuima at model@measure$intensity)))
idx.intensity <- append(idx.intensity,i)
}
-
+ }
+
assign("idx.intensity", idx.intensity, envir=env)
assign("measure.var", args[1], envir=env)
}
-
+
f <- function(p) {
mycoef <- as.list(p)
if(!is.CARMA(yuima)){
@@ -551,11 +553,11 @@
mycoef[fixed.par] <- fixed
minusquasilogl(yuima=yuima, param=mycoef, print=print, env,rcpp=rcpp)
}
-
+
# SMI-2/9/14:
fpsi <- function(p){
mycoef <- as.list(p)
-
+
idx.cont <- c(idx.diff,idx.drift)
if(length(c(idx.fixed,idx.cont))>0)
names(mycoef) <- nm[-c(idx.fixed,idx.cont)]
@@ -566,8 +568,8 @@
#print(p)
minusquasipsi(yuima=yuima, param=mycoef, print=print, env=env)
}
-
-
+
+
fj <- function(p) {
mycoef <- as.list(p)
# names(mycoef) <- nm
@@ -584,18 +586,18 @@
mycoef[fixed.par] <- fixed
minusquasilogl(yuima=yuima, param=mycoef, print=print, env,rcpp=rcpp)
}
-
+
oout <- NULL
HESS <- matrix(0, length(nm), length(nm))
colnames(HESS) <- nm
rownames(HESS) <- nm
-
-
+
+
HaveDriftHess <- FALSE
HaveDiffHess <- FALSE
HaveMeasHess <- FALSE
-
-
+
+
if(length(start)){
if(JointOptim){ ### joint optimization
old.fixed <- fixed
@@ -605,22 +607,22 @@
if(length(c(idx.fixed,idx.measure))>0)
new.start <- start[-c(idx.fixed,idx.measure)] # considering only initial guess for
}
-
+
if(length(new.start)>1){ #??multidimensional optim # Adjust lower for no negative Noise
if(is.CARMA(yuima) && (NoNeg.Noise==TRUE))
if(mean.noise %in% names(lower)){lower[mean.noise]<-10^-7}
oout <- optim(new.start, fj, method = method, hessian = TRUE, lower=lower, upper=upper)
-
+
if(length(fixed)>0)
oout$par[fixed.par]<- unlist(fixed)[fixed.par]
-
+
if(is.CARMA(yuima)){
HESS <- oout$hessian
} else {
HESS[names(new.start),names(new.start)] <- oout$hessian
}
-
-
+
+
if(is.CARMA(yuima) && length(yuima at model@info at scale.par)!=0){
b0<-paste0(yuima at model@info at ma.par,"0",collapse="")
idx.b0<-match(b0,rownames(HESS))
@@ -646,7 +648,7 @@
HESS<-HESS[,-indx.fixed]
}
}
-
+
if(is.CARMA(yuima) && length(yuima at model@parameter at measure)!=0){
for(i in c(1:length(fixed.par))){
indx.fixed<-match(fixed.par[i],rownames(HESS))
@@ -659,8 +661,8 @@
HESS<-HESS[,-idx.noise]
}
}
-
-
+
+
HaveDriftHess <- TRUE
HaveDiffHess <- TRUE
} else { ### one dimensional optim
@@ -674,13 +676,13 @@
} ### endif( length(start)>1 )
theta1 <- oout$par[diff.par]
theta2 <- oout$par[drift.par]
-
+
} else { ### first diffusion, then drift
theta1 <- NULL
-
+
old.fixed <- fixed
old.start <- start
-
+
if(length(idx.diff)>0){
## DIFFUSION ESTIMATIOn first
old.fixed <- fixed
@@ -694,7 +696,7 @@
fixed.par <- names(fixed)
idx.fixed <- match(fixed.par, nm)
names(new.start) <- nm[idx.diff]
-
+
mydots <- as.list(call)[-(1:2)]
mydots$print <- NULL
mydots$rcpp <- NULL #KK 08/07/16
@@ -708,7 +710,7 @@
mydots$joint <- NULL # LM 08/03/16
mydots$aggregation <- NULL # LM 08/03/16
mydots$threshold <- NULL #SMI 2/9/14
-
+
if((length(mydots$par)>1) | any(is.infinite(c(mydots$upper,mydots$lower)))){
mydots$method<-method ##song
oout <- do.call(optim, args=mydots)
@@ -727,16 +729,16 @@
oout <- list(par = theta1, value = opt1$objective)
}
theta1 <- oout$par
-
+
fixed <- old.fixed
start <- old.start
fixed.par <- old.fixed.par
-
+
} ## endif(length(idx.diff)>0)
-
+
theta2 <- NULL
-
-
+
+
if(length(idx.drift)>0){
## DRIFT estimation with first state diffusion estimates
fixed <- old.fixed
@@ -748,9 +750,9 @@
fixed <- new.fixed
fixed.par <- names(fixed)
idx.fixed <- match(fixed.par, nm)
-
+
names(new.start) <- nm[idx.drift]
-
+
mydots <- as.list(call)[-(1:2)]
mydots$print <- NULL
mydots$rcpp <- NULL #KK 08/07/16
@@ -757,7 +759,7 @@
mydots$fixed <- NULL
mydots$fn <- as.name("f")
mydots$threshold <- NULL #SMI 2/9/14
-
+
mydots$start <- NULL
mydots$par <- unlist(new.start)
mydots$hessian <- FALSE
@@ -765,9 +767,9 @@
mydots$lower <- unlist( lower[ nm[idx.drift] ])
mydots$joint <- NULL # LM 08/03/16
mydots$aggregation <- NULL # LM 08/03/16# LM 08/03/16
-
-
-
+
+
+
if(length(mydots$par)>1 | any(is.infinite(c(mydots$upper,mydots$lower)))){
if(is.CARMA(yuima)){
if(NoNeg.Noise==TRUE){
@@ -774,7 +776,7 @@
if((yuima at model@info at q+1)==yuima at model@info at p){
mydots$lower[names(start["NoNeg.Noise"])]<-10^(-7)
}
-
+
}
if(length(yuima at model@info at scale.par)!=0){
name_b0<-paste0(yuima at model@info at ma.par,"0",collapse="")
@@ -791,12 +793,12 @@
mydots$par <- unlist(new.start)
}
} # END if(is.CARMA)
-
+
mydots$method <- method #song
-
+
oout1 <- do.call(optim, args=mydots)
-
-
+
+
# oout1 <- optim(mydots$par,f,method = "L-BFGS-B" , lower = mydots$lower, upper = mydots$upper)
} else {
mydots$f <- mydots$fn
@@ -805,7 +807,7 @@
mydots$hessian <- NULL
mydots$method<-NULL
mydots$interval <- as.numeric(c(lower[drift.par],upper[drift.par]))
-
+
opt1 <- do.call(optimize, args=mydots)
theta2 <- opt1$minimum
names(theta2) <- drift.par
@@ -816,23 +818,23 @@
start <- old.start
old.fixed.par <- fixed.par
} ## endif(length(idx.drift)>0)
-
-
+
+
oout1 <- list(par= c(theta1, theta2))
if (! is.CARMA(yuima)){
if(length(c(diff.par, diff.par))>0)
names(oout1$par) <- c(diff.par,drift.par)
}
-
-
+
+
oout <- oout1
-
+
} ### endif JointOptim
} else {
list(par = numeric(0L), value = f(start))
}
-
-
+
+
fMeas <- function(p) {
mycoef <- as.list(p)
# if(! is.CARMA(yuima)){
@@ -842,8 +844,8 @@
minusquasipsi(yuima=yuima, param=mycoef, print=print, env=env)
# minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
}
-
-
+
+
fDrift <- function(p) {
mycoef <- as.list(p)
if(! is.CARMA(yuima)){
@@ -852,7 +854,7 @@
}
minusquasilogl(yuima=yuima, param=mycoef, print=print, env,rcpp=rcpp)
}
-
+
fDiff <- function(p) {
mycoef <- as.list(p)
if(! is.CARMA(yuima)){
@@ -861,37 +863,37 @@
}
minusquasilogl(yuima=yuima, param=mycoef, print=print, env,rcpp=rcpp)
}
-
+
# coef <- oout$par
#control=list()
#par <- coef
-
+
#names(par) <- unique(c(diff.par, drift.par))
# nm <- unique(c(diff.par, drift.par))
-
+
# START: ESTIMATION OF CP part
theta3 <- NULL
-
+
if(length(idx.measure)>0 & !is.CARMA(yuima)){
idx.cont <- c(idx.drift,idx.diff)
-
+
fixed <- old.fixed
start <- old.start
old.fixed.par <- fixed.par
new.fixed <- fixed
-
+
new.start <- start[idx.measure] # considering only initial guess for measure
new.fixed <- fixed
-
+
new.fixed[names(theta1)] <- theta1
new.fixed[names(theta2)] <- theta2
-
+
fixed <- new.fixed
fixed.par <- names(fixed)
idx.fixed <- match(fixed.par, nm)
# names(new.start) <- nm[idx.drift]
names(new.start) <- nm[idx.measure]
-
+
mydots <- as.list(call)[-(1:2)]
# mydots$print <- NULL
mydots$threshold <- NULL
@@ -899,7 +901,7 @@
mydots$fn <- as.name("fpsi")
mydots$start <- NULL
mydots$threshold <- NULL #SMI 2/9/14
-
+
mydots$par <- unlist(new.start)
mydots$hessian <- TRUE
mydots$joint <- NULL
@@ -906,36 +908,36 @@
mydots$upper <- unlist( upper[ nm[idx.measure] ])
mydots$lower <- unlist( lower[ nm[idx.measure] ])
mydots$method <- method
-
+
oout3 <- do.call(optim, args=mydots)
-
+
theta3 <- oout3$par
#print(theta3)
HESS[measure.par,measure.par] <- oout3$hessian
HaveMeasHess <- TRUE
-
+
fixed <- old.fixed
start <- old.start
fixed.par <- old.fixed.par
}
# END: ESTIMATION OF CP part
-
-
-
+
+
+
if(!is.CARMA(yuima)){
-
+
oout4 <- list(par= c(theta1, theta2, theta3))
names(oout4$par) <- c(diff.par,drift.par,measure.par)
oout <- oout4
}
-
+
coef <- oout$par
-
-
+
+
control=list()
par <- coef
if(!is.CARMA(yuima)){
-
+
names(par) <- unique(c(diff.par, drift.par,measure.par))
nm <- unique(c(diff.par, drift.par,measure.par))
} else {
@@ -943,19 +945,19 @@
nm <- unique(c(diff.par, drift.par))
}
#return(oout)
-
-
+
+
if(is.CARMA(yuima) && length(yuima at model@parameter at measure)!=0){
nm <-c(nm,measure.par)
if((NoNeg.Noise==TRUE)){nm <-c(nm,mean.noise)}
-
+
nm<-unique(nm)
}
if(is.CARMA(yuima) && (length(yuima at model@info at loc.par)!=0)){
nm <-unique(c(nm,yuima at model@info at loc.par))
}
-
-
+
+
conDrift <- list(trace = 5, fnscale = 1,
parscale = rep.int(5, length(drift.par)),
ndeps = rep.int(0.001, length(drift.par)), maxit = 100L,
@@ -989,9 +991,9 @@
beta = 0.5, gamma = 2, REPORT = 10, type = 1, lmm = 5,
factr = 1e+07, pgtol = 0, tmax = 10, temp = 10)
}
-
-
-
+
+
+
if(!HaveDriftHess & (length(drift.par)>0)){
#hess2 <- .Internal(optimhess(coef[drift.par], fDrift, NULL, conDrift))
if(!is.CARMA(yuima)){
@@ -1009,24 +1011,24 @@
HESS<-HESS[,-idx.b0]
}
}
-
+
if(!HaveDiffHess & (length(diff.par)>0)){
hess1 <- optimHess(coef[diff.par], fDiff, NULL, control=conDiff)
HESS[diff.par,diff.par] <- hess1
}
-
+
oout$hessian <- HESS
-
-
- if(!HaveMeasHess & (length(measure.par)>0) & !is.CARMA(yuima)){
+
+
+ if(!HaveMeasHess & (length(measure.par) > 0) & !is.CARMA(yuima)){
hess1 <- optimHess(coef[measure.par], fMeas, NULL, control=conMeas)
oout$hessian[measure.par,measure.par] <- hess1
}
-
+
# vcov <- if (length(coef))
# solve(oout$hessian)
# else matrix(numeric(0L), 0L, 0L)
-
+
vcov <- matrix(NA, length(coef), length(coef))
if (length(coef)) {
rrr <- try(solve(oout$hessian), TRUE)
@@ -1033,25 +1035,25 @@
if(class(rrr)[1] != "try-error")
vcov <- rrr
}
-
+
mycoef <- as.list(coef)
-
+
if(!is.CARMA(yuima)){
names(mycoef) <- nm
}
idx.fixed <- orig.idx.fixed
-
-
-
+
+
+
mycoef.cont <- mycoef
if(length(c(idx.fixed,idx.measure)>0)) # SMI 2/9/14
mycoef.cont <- mycoef[-c(idx.fixed,idx.measure)] # SMI 2/9/14
-
-
+
+
min.diff <- 0
min.jump <- 0
-
-
+
+
if(length(c(diff.par,drift.par))>0 & !is.CARMA(yuima)){ # LM 04/09/14
min.diff <- minusquasilogl(yuima=yuima, param=mycoef[c(diff.par,drift.par)], print=print, env,rcpp=rcpp)
}else{
@@ -1059,30 +1061,30 @@
min.diff <- minusquasilogl(yuima=yuima, param=mycoef, print=print, env,rcpp=rcpp)
}
}
-
+
if(length(c(measure.par))>0 & !is.CARMA(yuima))
min.jump <- minusquasipsi(yuima=yuima, param=mycoef[measure.par], print=print, env=env)
-
-
-
+
+
+
min <- min.diff + min.jump
if(min==0)
min <- NA
-
-
+
+
dummycov<-matrix(0,length(coef),length(coef))
rownames(dummycov)<-names(coef)
colnames(dummycov)<-names(coef)
dummycov[rownames(vcov),colnames(vcov)]<-vcov
vcov<-dummycov
-
-
+
+
# new("mle", call = call, coef = coef, fullcoef = unlist(mycoef),
# vcov = vcov, min = min, details = oout, minuslogl = minusquasilogl,
# method = method)
#LM 11/01
if(!is.CARMA(yuima)){
- if(length(measure.par)>0){
+ if(length(yuima at model@measure.type) > 0 && yuima at model@measure.type == "CP"){
final_res<-new("yuima.CP.qmle",
Jump.times=env$time[env$Cn.r==0],
Jump.values=env$deltaX[env$Cn.r==0,],
@@ -1116,35 +1118,35 @@
}
}
}
-
+
if(!is.CARMA(yuima)){
return(final_res)
}else {
-
+
param<-coef(final_res)
-
+
observ<-yuima at data
model<-yuima at model
info<-model at info
-
+
numb.ar<-info at p
name.ar<-paste(info at ar.par,c(numb.ar:1),sep="")
ar.par<-param[name.ar]
-
+
numb.ma<-info at q
name.ma<-paste(info at ma.par,c(0:numb.ma),sep="")
ma.par<-param[name.ma]
-
+
loc.par=NULL
if (length(info at loc.par)!=0){
loc.par<-param[info at loc.par]
}
-
+
scale.par=NULL
if (length(info at scale.par)!=0){
scale.par<-param[info at scale.par]
}
-
+
lin.par=NULL
if (length(info at lin.par)!=0){
lin.par<-param[info at lin.par]
@@ -1155,12 +1157,12 @@
yuima.warn("Insert constraints in Autoregressive parameters for enforcing stationarity" )
cat("\n Starting Estimation Increments ...\n")
}
-
+
ttt<-observ at zoo.data[[1]]
tt<-index(ttt)
y<-coredata(ttt)
if(NoNeg.Noise==TRUE && (info at p==(info at q+1))){final_res at coef[mean.noise]<-mean(y)/tail(ma.par,n=1)*ar.par[1]}
-
+
levy<-yuima.CarmaNoise(y,tt,ar.par,ma.par, loc.par, scale.par, lin.par, NoNeg.Noise)
inc.levy<-NULL
if (!is.null(levy)){
@@ -1176,23 +1178,23 @@
model = yuima at model, nobs=yuima.nobs, logL.Incr = NULL)
return(carma_final_res)
}
-
+
cat("\nStarting Estimation parameter Noise ...\n")
-
+
dummycovCarmapar<-vcov[unique(c(drift.par,diff.par)),unique(c(drift.par,diff.par))]
if(!is.null(loc.par)){
dummycovCarmapar<-vcov[unique(c(drift.par,diff.par,info at loc.par)),
unique(c(drift.par,diff.par,info at loc.par))]
}
-
-
-
+
+
+
dummycovCarmaNoise<-vcov[unique(measure.par),unique(c(measure.par))] #we need to adjusted
dummycoeffCarmapar<-coef[unique(c(drift.par,diff.par))]
if(!is.null(loc.par)){
dummycoeffCarmapar<-coef[unique(c(drift.par,diff.par,info at loc.par))]
}
-
+
dummycoeffCarmaNoise<-coef[unique(c(measure.par))]
coef<-NULL
coef<-c(dummycoeffCarmapar,dummycoeffCarmaNoise)
@@ -1200,7 +1202,7 @@
if(!is.null(loc.par)){
names.par<-c(unique(c(drift.par,diff.par,info at loc.par)),unique(c(measure.par)))
}
-
+
names(coef)<-names.par
cov<-NULL
cov<-matrix(0,length(names.par),length(names.par))
@@ -1211,9 +1213,9 @@
}else{
cov[unique(c(drift.par,diff.par,info at loc.par)),unique(c(drift.par,diff.par,info at loc.par))]<-dummycovCarmapar
}
-
+
cov[unique(c(measure.par)),unique(c(measure.par))]<-dummycovCarmaNoise
-
+
if(length(model at measure.type)!=0){
if(model at measure.type=="CP"){
name.func.dummy <- as.character(model at measure$df$expr[1])
@@ -1223,7 +1225,7 @@
name.int.dummy <- as.character(model at measure$intensity)
valueintensity<-as.numeric(name.int.dummy)
NaIdx<-which(!is.na(c(valueintensity,valuemeasure)))
-
+
if(length(NaIdx)!=0){
yuima.warn("the constrained MLE for levy increment will be implemented as soon as possible")
carma_final_res<-new("yuima.carma.qmle", call = call, coef = coef, fullcoef = unlist(mycoef),
@@ -1232,7 +1234,7 @@
model = yuima at model, logL.Incr = NULL)
return(carma_final_res)
}
-
+
if(aggregation==TRUE){
if(floor(yuima at sampling@n/yuima at sampling@Terminal)!=yuima at sampling@n/yuima at sampling@Terminal){
yuima.stop("the n/Terminal in sampling information is not an integer. Set Aggregation=FALSE")
@@ -1244,16 +1246,16 @@
}else{
inc.levy1<-inc.levy
}
-
+
names.measpar<-c(name.int.dummy, names.measpar)
-
+
if(measurefunc=="dnorm"){
-
+
# result.Lev<-yuima.Estimation.CPN(Increment.lev=inc.levy1,param0=coef[ names.measpar],
# fixed.carma=fixed.carma,
# lower.carma=lower.carma,
# upper.carma=upper.carma)
-
+
result.Lev<-yuima.Estimation.Lev(Increment.lev=inc.levy1,
param0=coef[ names.measpar],
fixed.carma=fixed.carma,
@@ -1263,7 +1265,7 @@
measure.type=model at measure.type,
dt=env$h,
aggregation=aggregation)
-
+
}
if(measurefunc=="dgamma"){
# result.Lev<-yuima.Estimation.CPGam(Increment.lev=inc.levy1,param0=coef[ names.measpar],
@@ -1270,7 +1272,7 @@
# fixed.carma=fixed.carma,
# lower.carma=lower.carma,
# upper.carma=upper.carma)
-
+
result.Lev<-yuima.Estimation.Lev(Increment.lev=inc.levy1,
param0=coef[ names.measpar],
fixed.carma=fixed.carma,
@@ -1286,7 +1288,7 @@
# fixed.carma=fixed.carma,
# lower.carma=lower.carma,
# upper.carma=upper.carma)
-
+
result.Lev<-yuima.Estimation.Lev(Increment.lev=inc.levy1,
param0=coef[ names.measpar],
fixed.carma=fixed.carma,
@@ -1296,18 +1298,18 @@
measure.type=model at measure.type,
dt=env$h,
aggregation=aggregation)
-
+
}
Inc.Parm<-result.Lev$estLevpar
IncVCOV<-result.Lev$covLev
-
+
names(Inc.Parm)[NaIdx]<-measure.par
rownames(IncVCOV)[NaIdx]<-as.character(measure.par)
colnames(IncVCOV)[NaIdx]<-as.character(measure.par)
-
+
coef<-NULL
coef<-c(dummycoeffCarmapar,Inc.Parm)
-
+
names.par<-names(coef)
cov<-NULL
cov<-matrix(0,length(names.par),length(names.par))
@@ -1319,8 +1321,8 @@
cov[unique(c(drift.par,diff.par,info at loc.par)),unique(c(drift.par,diff.par,info at loc.par))]<-dummycovCarmapar
}
cov[names(Inc.Parm),names(Inc.Parm)]<-IncVCOV
-
-
+
+
}
if(yuima at model@measure.type=="code"){
# # "rIG", "rNIG", "rgamma", "rbgamma", "rvgamma"
@@ -1332,7 +1334,7 @@
name.func<- substr(name.func.dummy,1,(nchar(name.func.dummy)-1))
names.measpar<-as.vector(strsplit(name.func,', '))[[1]][-1]
valuemeasure<-as.numeric(names.measpar)
-
+
NaIdx<-which(!is.na(valuemeasure))
}
if(length(NaIdx)!=0){
@@ -1355,7 +1357,7 @@
inc.levy1<-inc.levy
}
if(measurefunc=="yuima.law"){
-
+
dummyParMeas<-c(coef[measure.par],1)
names(dummyParMeas)<-c(measure.par,yuima at model@time.variable)
cond <- length(dens(yuima at model@measure$df,x=as.numeric(inc.levy1),param=as.list(dummyParMeas)))
@@ -1364,7 +1366,7 @@
covLev=matrix(NA,
length(coef[measure.par]),
length(coef[measure.par]))
- )
+ )
yuima.warn("Levy measure parameters can not be estimated.")
}else{
dummyMyfunMeas<-function(par, Law, Data, time, param.name, name.time){
@@ -1388,21 +1390,21 @@
}
prova <- optim(fn = dummyMyfunMeas, par = coef[measure.par],
- method = method,Law=yuima at model@measure$df,
- Data=inc.levy1,
- time=mytime, param.name=measure.par,
- name.time = yuima at model@time.variable)
+ method = method,Law=yuima at model@measure$df,
+ Data=inc.levy1,
+ time=mytime, param.name=measure.par,
+ name.time = yuima at model@time.variable)
Heeee<-optimHess(fn = dummyMyfunMeas, par = coef[measure.par],
- Law=yuima at model@measure$df,
- Data=inc.levy1,
- time=mytime, param.name=measure.par,
- name.time = yuima at model@time.variable)
+ Law=yuima at model@measure$df,
+ Data=inc.levy1,
+ time=mytime, param.name=measure.par,
+ name.time = yuima at model@time.variable)
result.Lev <- list(estLevpar=prova$par,covLev=solve(Heeee))
}
}
-
+
if(measurefunc=="rIG"){
-
+
# result.Lev<-list(estLevpar=coef[ names.measpar],
# covLev=matrix(NA,
# length(coef[ names.measpar]),
@@ -1412,7 +1414,7 @@
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 751
More information about the Yuima-commits
mailing list