[Yuima-commits] r394 - in pkg/yuima: . R man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Nov 19 21:21:23 CET 2015
Author: lorenzo
Date: 2015-11-19 21:21:23 +0100 (Thu, 19 Nov 2015)
New Revision: 394
Added:
pkg/yuima/R/PseudoLogLikCOGARCH.R
pkg/yuima/src/PseudoLoglikCOGARCH.c
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/R/qmle.R
pkg/yuima/R/simulate.R
pkg/yuima/man/qmle.Rd
Log:
Added methods and Class for estimation of a COGARCH(p,q) model using the maximization of the pseudo-likelihood function
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2015-10-10 06:36:25 UTC (rev 393)
+++ pkg/yuima/DESCRIPTION 2015-11-19 20:21:23 UTC (rev 394)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project Package for SDEs
-Version: 1.0.75
+Version: 1.0.76
Depends: R(>= 2.10.0), methods, zoo, stats4, utils, expm, cubature, mvtnorm
Author: YUIMA Project Team
Maintainer: Stefano M. Iacus <stefano.iacus at unimi.it>
Added: pkg/yuima/R/PseudoLogLikCOGARCH.R
===================================================================
--- pkg/yuima/R/PseudoLogLikCOGARCH.R (rev 0)
+++ pkg/yuima/R/PseudoLogLikCOGARCH.R 2015-11-19 20:21:23 UTC (rev 394)
@@ -0,0 +1,316 @@
+# This code is useful for estimation of the COGARCH(p,q) model according the PseudoLogLikelihood
+# maximization procedure developed in Iacus et al. 2015
+#
+
+PseudoLogLik.COGARCH <- function(yuima, start, method="BFGS", fixed = list(),
+ lower, upper, Est.Incr, call, grideq, ...){
+
+ if(is(yuima,"yuima")){
+ model <- yuima at model
+ info <- model at info
+ Data <- onezoo(yuima)
+ }
+
+ time <- index(Data)
+ Obs <- as.numeric(as.matrix(Data)[,1])
+ my.env <- new.env()
+ param <- unlist(start)
+
+ meas.par <- model at parameter@measure
+
+ if(length(meas.par)==0 && Est.Incr=="IncrPar"){
+ yuima.warn("The dimension of measure parameters is zero, yuima changes 'Est.Incr = IncrPar' into 'Est.Incr = Incr'")
+ Est.Incr <- "Incr"
+ }
+
+ ar.names <- paste0(info at ar.par,c(1:info at q))
+ ma.names <- paste0(info at ma.par,c(1:info at p))
+
+ start.state <- paste0(paste0(info at Latent.var,0),c(1:info at q))
+
+ e <- matrix(0, info at q,1)
+ e[info at q,1] <- 1
+ assign("e", e, envir = my.env)
+
+ assign("start.state", start.state, envir = my.env)
+ assign("q", info at q, envir = my.env)
+ assign("p", info at p, envir = my.env)
+ assign("B", matrix(0,info at q,info at q), envir = my.env)
+ # consider two cases:
+ # 1) equally spaced grid
+ if(grideq){
+ assign("Deltat", tail(index(Data),1)/length(index(Data)), envir = my.env)
+ # assign("Deltat", diff(time)[1], envir = my.env)
+ # 2) no-equally spaced grid
+ assign("grideq", TRUE, envir = my.env)
+ }else{
+ assign("Deltat", diff(time), envir = my.env)
+ assign("grideq", FALSE, envir = my.env)
+ }
+ assign("Obs", (diff(Obs))^2, envir = my.env)
+ assign("nObs",length(Obs),envir = my.env)
+ assign("ar.names", ar.names, envir = my.env)
+ assign("ma.names", ma.names, envir = my.env)
+ assign("loc.par",info at loc.par, envir = my.env)
+
+ I <- diag(info at q)
+ assign("I",I, envir = my.env)
+
+ out<-NULL
+ if(length(lower)==0 && length(upper)>0 && length(fixed)==0){
+ out <- optim(par=param, fn=minusloglik.COGARCH1,
+ method = method, upper=upper, env = my.env,...)
+
+
+ }
+
+ if(length(lower)==0 && length(upper)==0 && length(fixed)>0){
+ out <- optim(par=param, fn=minusloglik.COGARCH1,
+ method = method, fixed=fixed, env = my.env,...)
+
+ }
+
+
+ if(length(lower)>0 && length(upper)==0 && length(fixed)==0){
+ out <- optim(par=param, fn=minusloglik.COGARCH1,
+ method = method, lower=lower, env = my.env,...)
+ }
+
+ if(length(lower)>0 && length(upper)>0 && length(fixed)==0){
+ out <- optim(par=param, fn=minusloglik.COGARCH1,
+ method = method, upper = upper,
+ lower=lower, env = my.env,...)
+ }
+
+
+ if(length(lower)==0 && length(upper)>0 && length(fixed)>0){
+ out <- optim(par=param, fn=minusloglik.COGARCH1,
+ method = method, upper = upper,
+ fixed = fixed, env = my.env,...)
+ }
+
+ if(length(lower)>0 && length(upper)==0 && length(fixed)>0){
+ out <- optim(par=param, fn=minusloglik.COGARCH1,
+ method = method, lower = lower,
+ fixed = fixed, env = my.env,...)
+ }
+
+
+ if(length(lower)>0 && length(upper)>0 && length(fixed)>0){
+ out <- optim(par=param, fn=minusloglik.COGARCH1,
+ method = method, lower = lower,
+ fixed = fixed, upper = upper,
+ env = my.env,...)
+ }
+
+
+ if(is.null(out)){
+ out <- optim(par=param, fn=minusloglik.COGARCH1,
+ method = method, env = my.env,...)
+ }
+
+# control= list(maxit=100),
+ # Write the object mle with result
+
+
+ bvect<-out$par[ar.names]
+ bq<-bvect[1]
+ avect<-out$par[ma.names]
+ a1<-avect[1]
+
+ a0<-out$par[info at loc.par]
+
+ if(length(meas.par)!=0){
+ idx.dumm<-match(meas.par,names(out$par))
+ out$par<-out$par[- idx.dumm]
+ }
+
+ idx.dumm1<-match(start.state,names(out$par))
+ coef <- out$par[-idx.dumm1]
+
+ vcov<-matrix(NA, length(coef), length(coef))
+ names_coef<-names(coef)
+ colnames(vcov) <- names_coef
+ rownames(vcov) <- names_coef
+ mycoef <- start
+ # min <- out$value
+ objFun <- "PseudoLogLik"
+ min <- numeric()
+
+
+ res<-new("cogarch.gmm", call = call, coef = coef, fullcoef = unlist(coef),
+ vcov = vcov, min = min, details = list(),
+ method = character(),
+ model = model,
+ objFun = objFun
+ )
+
+
+ return(res)
+}
+
+minusloglik.COGARCH1<-function(param,env){
+
+# assign("start.state", start.state, envir = my.env)
+# assign("q", info at q, envir = my.env)
+# assign("p", info at p, envir = my.env)
+# assign("time", time, envir = my.env)
+# assign("Obs", Obs, envir = my.env)
+# assign("nObs",length(Obs),envir = my.env)
+# assign("ar.names", ar.names, envir = my.env)
+# assign("ma.names", ma.names, envir = my.env)
+# assign("loc.par",info at loc.par, envir = my.env)
+ a0<-param[env$loc.par]
+ bq<-param[env$ar.names[env$q]]
+ a1<-param[env$ma.names[1]]
+
+ stateMean <- a0/(bq-a1)*as.matrix(c(1,numeric(length=(env$q-1))))
+
+ param[env$start.state]<-stateMean
+ state <- stateMean
+# state <- param[env$start.state]
+ DeltaG2 <- env$Obs
+ B <- env$B
+ if(env$q>1){
+ B[1:(env$q-1),] <- c(matrix(0,(env$p-1),1), diag(env$q-1))
+ }
+ B[env$q,] <- -param[env$ar.names[env$q:1]]
+ a<-matrix(0,env$q,1)
+ a[1:env$p,]<-param[env$ma.names]
+
+ ta<-t(a)
+ e <- env$e
+ Btilde <-B+e%*%ta
+ InvBtilde <- solve(Btilde)
+ V1 <- a0+ta%*% state
+ V <- V1[1]
+ Deltat<- env$Deltat
+ I <- env$I
+# VarDeltaG <- 0
+ PseudologLik <- 0
+
+# DeltatB1 <- lapply(as.list(Deltat), function(dt,B){expm(B*dt)} , B)
+# # DeltatB <- lapply(as.list(Deltat), "*" , B)
+# # assign("DeltatB",as.list(DeltatB),.GlobalEnv)
+# outputB <- matrix(unlist(DeltatB1), ncol = env$q, byrow = TRUE)
+ if(env$grideq){
+ DeltatB1 <- expm(B*Deltat)
+
+# DeltatB2 <- lapply(as.list(Deltat), function(dt,B){expm(B*dt)} , Btilde)
+# # DeltatB <- lapply(as.list(Deltat), "*" , B)
+# # assign("DeltatB",as.list(DeltatB),.GlobalEnv)
+# outputB2 <- matrix(unlist(DeltatB2), ncol = env$q, byrow = TRUE)
+
+ DeltatB2 <- expm(Btilde*Deltat)
+
+# DeltatB3 <- lapply(as.list(-Deltat), function(dt,B){expm(B*dt)} , Btilde)
+# outputB3 <- matrix(unlist(DeltatB3), ncol = env$q, byrow = TRUE)
+
+ DeltatB3 <- expm(-Btilde*Deltat)
+
+ dummyMatr <- ta%*%DeltatB2%*%InvBtilde%*%(I-DeltatB3)
+ dummyeB1 <- e%*%ta%*%DeltatB1
+
+ #aa <- .Call("myfun1", DeltatB1, state)
+ PseudologLik <-.Call("pseudoLoglik_COGARCH1", a0, bq, a1, stateMean, Q=as.integer(env$q),
+ DeltaG2, Deltat, DeltatB1, a, e,
+ V, nObs=as.integer(env$nObs-1),
+ dummyMatr, dummyeB1)
+ #cat(sprintf("\n%.5f ", PseudologLik))
+#
+#
+# PseudologLikR<-0
+# V1 <- a0+ta%*% state
+# V <- V1[1]
+#
+#
+# for(i in c(1:(env$nObs-1))){
+#
+# # cat(sprintf("\n dummy1R %.10f ",dummyMatr%*%(state-stateMean) ))
+# # d <- 0
+# # for(j in 1:2){
+# # d<- d+dummyMatr[1,j]*(state[j]-stateMean[j])
+# # cat(sprintf("\n %d dummy1R %.10f ",j,d ))
+# # }
+# VarDeltaG <- a0*Deltat*bq/(bq-a1)+ dummyMatr%*%(state-stateMean)
+#
+# VarDeltaG <- VarDeltaG[1]
+# #state <- (I+DeltaG2[i]/V*e%*%ta)%*%DeltatB1%*%state+a0*DeltaG2[i]/V*e
+# state <- DeltatB1%*%state+DeltaG2[i]/V*dummyeB1%*%state+a0*DeltaG2[i]/V*e
+# #state <- DeltatB1%*%state+dummyeB1%*%state
+# # cat(sprintf("\n d1 %.10f d2 %.10f", DeltatB1%*%state, dummyeB1%*%state));
+# V <- a0+ta%*% state
+# V <- V[1]
+# if(is.nan(VarDeltaG))
+# VarDeltaG<- 10^(-6)
+# PseudologLikR<--0.5*(DeltaG2[i]/VarDeltaG+log(VarDeltaG)+log(2.*3.14159265))+PseudologLikR
+# # cat(sprintf("\n%.5f - %.5f %.5f - %.5f",VarDeltaG, state[1], state[2],V ))
+# cat(sprintf("\n Part %.10f partial %.10f ", PseudologLikR, VarDeltaG))
+# if(is.nan(V))
+# V <- 10^(-6)
+#
+# }
+# #
+# cat(sprintf("\n%.5f - %.5f",PseudologLikR, PseudologLik))
+# # #
+# #
+
+ }else{
+ DeltatB1 <- lapply(as.list(Deltat), function(dt,B){expm(B*dt)} , B)
+ DeltatB2 <- lapply(as.list(Deltat), function(dt,B){expm(B*dt)} , Btilde)
+ DeltatB3 <- lapply(as.list(-Deltat), function(dt,B){expm(B*dt)} , Btilde)
+
+ for(i in c(1:(env$nObs-1))){
+ VarDeltaG <- as.numeric(a0*Deltat[i]*bq/(bq-a1)+ta%*%DeltatB2[[i]]%*%InvBtilde%*%(I-DeltatB3[[i]])%*%(state-stateMean))
+ state <- (I+DeltaG2[i]/V*e%*%ta)%*%DeltatB1[[i]]%*%state+a0*DeltaG2[i]/V*e
+ V <- as.numeric(a0+ta%*% state)
+ PseudologLik<--1/2*(DeltaG2[i]/VarDeltaG+log(VarDeltaG)+log(2*pi))+PseudologLik
+ }
+ }
+#
+# PseudologLik <- 0
+#
+# for(i in c(1:(env$nObs-1))){
+# VarDeltaG <- a0*Deltat*bq/(bq-a1)+ dummyMatr%*%(state-stateMean)
+# #state <- (I+DeltaG2[i]/V*e%*%ta)%*%DeltatB1%*%state+a0*DeltaG2[i]/V*e
+# state <- DeltatB1%*%state+DeltaG2[i]/V*dummyeB1%*%state+a0*DeltaG2[i]/V*e
+# V <- as.numeric(a0+ta%*% state)
+# PseudologLik<--1/2*(DeltaG2[i]/VarDeltaG+log(VarDeltaG)+log(2*pi))+PseudologLik
+# }
+
+# for(i in c(1:(env$nObs-1))){
+# VarDeltaG <- as.numeric(a0*Deltat[i]*bq/(bq-a1)+t(a)%*%DeltatB2[[i]]%*%InvBtilde%*%(I-DeltatB3[[i]])%*%(state-stateMean))
+
+
+
+ # state <- (I+DeltaG2[i]/V*e%*%t(a))%*%DeltatB1[[i]]%*%state+a0*DeltaG2[i]/V*e
+# V <- as.numeric(a0+t(a)%*% state)
+# PseudologLik<--1/2*(DeltaG2[i]/VarDeltaG+log(VarDeltaG)+log(2*pi))
+# }
+
+# dummyMatr <- ta%*%DeltatB2%*%InvBtilde%*%(I-DeltatB3)
+# dummyeB1 <- e%*%ta%*%DeltatB1
+# PseudologLik1 <- 0
+# for(i in c(1:(env$nObs-1))){
+# VarDeltaG <- as.numeric(a0*Deltat*bq/(bq-a1)+dummyMatr%*%(state-stateMean))
+# state <- DeltatB1%*%state+DeltaG2[i]/V*dummyeB1%*%state+a0*DeltaG2[i]/V*e
+# V <- as.numeric(a0+ta%*% state)
+# PseudologLik1 <- -1/2*(DeltaG2[i]/VarDeltaG+log(VarDeltaG)+log(2*pi))
+# if(is.finite(PseudologLik1)){
+# PseudologLik <- PseudologLik1 + PseudologLik
+# }
+# }
+
+ minusPseudoLogLik <- -PseudologLik
+ return(minusPseudoLogLik)
+}
+
+# res<-.Call("pseudoLoglik_COGARCH", a0, bq, a1, stateMean, Q=as.integer(env$q),
+# state, DeltaG2, Deltat, DeltatB, B, a, e,
+# Btilde, InvBtilde, V, I, VarDeltaG,
+# PseudologLik, nObs = as.integer(env$nObs-1), fn = quote(expm(x)) , rho= .GlobalEnv,
+# PACKAGE = "yuima")
+#
+# output <- matrix(unlist(res), ncol = env$q, byrow = TRUE)
+# res <- res
+
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2015-10-10 06:36:25 UTC (rev 393)
+++ pkg/yuima/R/qmle.R 2015-11-19 20:21:23 UTC (rev 394)
@@ -6,7 +6,7 @@
### TO BE FIXED: all caculations should be made on a private environment to
### avoid problems.
### I have rewritten drift.term and diff.term instead of calc.drift and
-### calc.diffusion to make them independent of the specification of the
+### calc.diffusion to make them independent of the specification of the
### parameters. S.M.I. 22/06/2010
drift.term <- function(yuima, theta, env){
@@ -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)
}
@@ -33,7 +33,7 @@
drift[,d] <- eval(DRIFT[d], envir=tmp.env)
}
- return(drift)
+ return(drift)
}
@@ -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)
}
@@ -131,43 +131,61 @@
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(),
+ lower, upper, Est.Incr, call, ...)
+ }else{
+ res <- PseudoLogLik.COGARCH(yuima, start, method=method, fixed = list(),
+ lower, upper, Est.Incr, call, grideq = FALSE,...)
+ }
+
+ 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) )
+## 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
+ # 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")
@@ -176,7 +194,7 @@
measurefunc <- substring(yuima at model@measure$df$exp, 1, tmp-1)
if(!is.na(match(measurefunc,CPlist))){
yuima.warn("carma(p,q): the qmle for a carma(p,q) driven by a Compound Poisson with no-negative random size")
- NoNeg.Noise<-TRUE
+ NoNeg.Noise<-TRUE
# we need to add a new parameter for the mean of the Noise
if((yuima at model@info at q+1)==(yuima at model@info at q+1)){
start[["mean.noise"]]<-1
@@ -185,7 +203,7 @@
}
}
-
+
if(yuima at model@measure.type=="code"){
tmp <- regexpr("\\(", yuima at model@measure$df$exp)[1]
measurefunc <- substring(yuima at model@measure$df$exp, 1, tmp-1)
@@ -196,31 +214,31 @@
start[["mean.noise"]]<-1
}
#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)
+ 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
@@ -243,9 +261,9 @@
JointOptim <- TRUE
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
@@ -262,25 +280,25 @@
# if(length(jump.par)+length(measure.par)>0)
# 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
+ # 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)){
@@ -289,14 +307,14 @@
fullcoef<-c(fullcoef, mean.noise)
}
}
-
+
# if(is.CARMA(yuima) && (length(measure.par)>0)){
fullcoef<-c(fullcoef, measure.par)
#}
npar <- length(fullcoef)
-
-
+
+
fixed.par <- names(fixed) # We use Fixed.par when we consider a Carma with scale parameter
if(is.CARMA(yuima) && (length(measure.par)>0)){
fixed.carma=NULL
@@ -332,30 +350,30 @@
}
}
}
-
-
-
+
+
+
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])
fixed[measure.par[j]]<-start[measure.par[j]]
}
}
-
+
}
- if (any(!(fixed.par %in% fullcoef)))
+ 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
@@ -366,50 +384,50 @@
idx.xinit <- as.integer(na.omit(match(xinit.par,nm)))# We need to add idx if NoNeg.Noise is TRUE
#}
}
-
+
idx.fixed <- match(fixed.par, nm)
orig.idx.fixed <- idx.fixed
-
+
tmplower <- as.list( rep( -Inf, length(nm)))
- names(tmplower) <- nm
+ names(tmplower) <- nm
if(!missing(lower)){
idx <- match(names(lower), names(tmplower))
if(any(is.na(idx)))
yuima.stop("names in 'lower' do not match names fo parameters")
- tmplower[ idx ] <- lower
+ tmplower[ idx ] <- lower
}
lower <- tmplower
-
+
tmpupper <- as.list( rep( Inf, length(nm)))
- names(tmpupper) <- nm
+ names(tmpupper) <- nm
if(!missing(upper)){
idx <- match(names(upper), names(tmpupper))
if(any(is.na(idx)))
yuima.stop("names in 'lower' do not match names fo parameters")
- tmpupper[ idx ] <- upper
+ tmpupper[ idx ] <- upper
}
upper <- tmpupper
-
-
-
-
+
+
+
+
d.size <- yuima at model@equation.number
if (is.CARMA(yuima)){
# 24/12
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)
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
@@ -423,33 +441,33 @@
# env$time.obs<-length(env$X)
#p <-yuima at model@info at p
assign("p", yuima at model@info at p, envir=env)
- assign("q", yuima at model@info at q, 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)
# p <-yuima at model@info at p
# 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)
-
+ 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)
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){
-
-
+
+
args <- unlist(strsplit(suppressWarnings(sub("^.+?\\((.+)\\)", "\\1",yuima at model@measure$df$expr,perl=TRUE)), ","))
idx.intensity <- numeric(0)
for(i in 1:length(measure.par)){
@@ -478,11 +496,11 @@
mycoef[fixed.par] <- fixed
minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
}
-
+
# 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)]
@@ -493,8 +511,8 @@
#print(p)
minusquasipsi(yuima=yuima, param=mycoef, print=print, env=env)
}
-
-
+
+
fj <- function(p) {
mycoef <- as.list(p)
# names(mycoef) <- nm
@@ -516,13 +534,13 @@
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
@@ -543,8 +561,8 @@
} 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))
@@ -575,9 +593,9 @@
} else { ### first diffusion, then drift
theta1 <- NULL
- old.fixed <- fixed
+ old.fixed <- fixed
old.start <- start
-
+
if(length(idx.diff)>0){
## DIFFUSION ESTIMATIOn first
old.fixed <- fixed
@@ -585,13 +603,13 @@
old.fixed.par <- fixed.par
new.start <- start[idx.diff] # considering only initial guess for diffusion
new.fixed <- fixed
- if(length(idx.drift)>0)
+ if(length(idx.drift)>0)
new.fixed[nm[idx.drift]] <- start[idx.drift]
fixed <- new.fixed
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$fixed <- NULL
@@ -602,15 +620,15 @@
mydots$upper <- as.numeric(unlist( upper[ nm[idx.diff] ]))
mydots$lower <- as.numeric(unlist( lower[ nm[idx.diff] ]))
mydots$threshold <- NULL #SMI 2/9/14
-
+
if((length(mydots$par)>1) | any(is.infinite(c(mydots$upper,mydots$lower)))){
oout <- do.call(optim, args=mydots)
} else {
mydots$f <- mydots$fn
mydots$fn <- NULL
mydots$par <- NULL
- mydots$hessian <- NULL
- mydots$method <- NULL
+ mydots$hessian <- NULL
+ mydots$method <- NULL
mydots$interval <- as.numeric(c(unlist(lower[diff.par]),unlist(upper[diff.par])))
@@ -619,19 +637,19 @@
opt1 <- do.call(optimize, args=mydots)
theta1 <- opt1$minimum
names(theta1) <- diff.par
- oout <- list(par = theta1, value = opt1$objective)
+ 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
@@ -643,7 +661,7 @@
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)]
@@ -651,24 +669,24 @@
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
mydots$upper <- unlist( upper[ nm[idx.drift] ])
mydots$lower <- unlist( lower[ nm[idx.drift] ])
-
-
-
-
-
+
+
+
+
+
if(length(mydots$par)>1 | any(is.infinite(c(mydots$upper,mydots$lower)))){
if(is.CARMA(yuima)){
if(NoNeg.Noise==TRUE){
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="")
@@ -685,47 +703,47 @@
mydots$par <- unlist(new.start)
}
} # END if(is.CARMA)
-
-
-
+
+
+
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
mydots$fn <- NULL
mydots$par <- NULL
- mydots$hessian <- NULL
- mydots$method <- NULL
- mydots$interval <- as.numeric(c(lower[drift.par],upper[drift.par]))
+ 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
- oout1 <- list(par = theta2, value = as.numeric(opt1$objective))
+ oout1 <- list(par = theta2, value = as.numeric(opt1$objective))
}
theta2 <- oout1$par
fixed <- old.fixed
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)){
@@ -754,7 +772,7 @@
}
minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
}
-
+
# coef <- oout$par
#control=list()
#par <- coef
@@ -764,7 +782,7 @@
# START: ESTIMATION OF CP part
theta3 <- NULL
-
+
if(length(idx.measure)>0 & !is.CARMA(yuima)){
idx.cont <- c(idx.drift,idx.diff)
@@ -775,16 +793,16 @@
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
@@ -792,43 +810,43 @@
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
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 {
@@ -836,12 +854,12 @@
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)){
@@ -850,16 +868,16 @@
conDrift <- list(trace = 5, fnscale = 1,
- parscale = rep.int(5, length(drift.par)),
- ndeps = rep.int(0.001, length(drift.par)), maxit = 100L,
- abstol = -Inf, reltol = sqrt(.Machine$double.eps), alpha = 1,
- beta = 0.5, gamma = 2, REPORT = 10, type = 1, lmm = 5,
+ parscale = rep.int(5, length(drift.par)),
+ ndeps = rep.int(0.001, length(drift.par)), maxit = 100L,
+ abstol = -Inf, reltol = sqrt(.Machine$double.eps), alpha = 1,
+ beta = 0.5, gamma = 2, REPORT = 10, type = 1, lmm = 5,
factr = 1e+07, pgtol = 0, tmax = 10, temp = 10)
- conDiff <- list(trace = 5, fnscale = 1,
- parscale = rep.int(5, length(diff.par)),
- ndeps = rep.int(0.001, length(diff.par)), maxit = 100L,
- abstol = -Inf, reltol = sqrt(.Machine$double.eps), alpha = 1,
- beta = 0.5, gamma = 2, REPORT = 10, type = 1, lmm = 5,
+ conDiff <- list(trace = 5, fnscale = 1,
+ parscale = rep.int(5, length(diff.par)),
+ ndeps = rep.int(0.001, length(diff.par)), maxit = 100L,
+ abstol = -Inf, reltol = sqrt(.Machine$double.eps), alpha = 1,
+ beta = 0.5, gamma = 2, REPORT = 10, type = 1, lmm = 5,
factr = 1e+07, pgtol = 0, tmax = 10, temp = 10)
conMeas <- list(trace = 5, fnscale = 1,
parscale = rep.int(5, length(measure.par)),
@@ -868,18 +886,18 @@
beta = 0.5, gamma = 2, REPORT = 10, type = 1, lmm = 5,
factr = 1e+07, pgtol = 0, tmax = 10, temp = 10)
if(is.CARMA(yuima) && length(yuima at model@info at loc.par)!=0 ){
- conDrift <- list(trace = 5, fnscale = 1,
- parscale = rep.int(5, length(c(drift.par,yuima at model@info at loc.par))),
- ndeps = rep.int(0.001, length(c(drift.par,yuima at model@info at loc.par))),
- maxit = 100L,
- abstol = -Inf, reltol = sqrt(.Machine$double.eps), alpha = 1,
- beta = 0.5, gamma = 2, REPORT = 10, type = 1, lmm = 5,
+ conDrift <- list(trace = 5, fnscale = 1,
+ parscale = rep.int(5, length(c(drift.par,yuima at model@info at loc.par))),
+ ndeps = rep.int(0.001, length(c(drift.par,yuima at model@info at loc.par))),
+ maxit = 100L,
+ abstol = -Inf, reltol = sqrt(.Machine$double.eps), alpha = 1,
+ beta = 0.5, gamma = 2, REPORT = 10, type = 1, lmm = 5,
factr = 1e+07, pgtol = 0, tmax = 10, temp = 10)
- conDiff <- list(trace = 5, fnscale = 1,
- parscale = rep.int(5, length(diff.par)),
- ndeps = rep.int(0.001, length(diff.par)), maxit = 100L,
- abstol = -Inf, reltol = sqrt(.Machine$double.eps), alpha = 1,
- beta = 0.5, gamma = 2, REPORT = 10, type = 1, lmm = 5,
+ conDiff <- list(trace = 5, fnscale = 1,
+ parscale = rep.int(5, length(diff.par)),
+ ndeps = rep.int(0.001, length(diff.par)), maxit = 100L,
+ abstol = -Inf, reltol = sqrt(.Machine$double.eps), alpha = 1,
+ beta = 0.5, gamma = 2, REPORT = 10, type = 1, lmm = 5,
factr = 1e+07, pgtol = 0, tmax = 10, temp = 10)
}
@@ -902,15 +920,15 @@
HESS<-HESS[,-idx.b0]
}
}
-
+
if(!HaveDiffHess & (length(diff.par)>0)){
hess1 <- optimHess(coef[diff.par], fDiff, NULL, control=conDiff)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 394
More information about the Yuima-commits
mailing list