[Yuima-commits] r267 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jan 13 22:07:49 CET 2014
Author: lorenzo
Date: 2014-01-13 22:07:48 +0100 (Mon, 13 Jan 2014)
New Revision: 267
Added:
pkg/yuima/R/CarmaRecovNoise.R
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/R/qmle.R
pkg/yuima/man/qmle.Rd
Log:
added invisible function CarmaRecovNoise.R for recovering the underlying noise for carma model.
Modified qmle function for Carma driven by a Levy
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2014-01-12 16:34:58 UTC (rev 266)
+++ pkg/yuima/DESCRIPTION 2014-01-13 21:07:48 UTC (rev 267)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package (unstable version)
-Version: 0.1.219
-Date: 2013-12-26
+Version: 0.1.220
+Date: 2014-01-13
Depends: methods, zoo, stats4, utils, Matrix
Suggests: cubature, mvtnorm
Author: YUIMA Project Team.
Added: pkg/yuima/R/CarmaRecovNoise.R
===================================================================
--- pkg/yuima/R/CarmaRecovNoise.R (rev 0)
+++ pkg/yuima/R/CarmaRecovNoise.R 2014-01-13 21:07:48 UTC (rev 267)
@@ -0,0 +1,231 @@
+
+# In this file we develop the procedure described in Brockwell, Davis and Yang (2012)
+# for estimation of the underlying noise once the parameters of carma(p,q) are previously
+# obtained
+
+
+yuima.carma.eigen<-function(A){
+ diagA<-eigen(A)
+ diagA$values<-diagA$values[order(diagA$values, na.last = TRUE, decreasing = TRUE)]
+ n_eigenval<-length(diagA$values)
+ diagA$vectors<-matrix(diagA$values[1]^(c(1:n_eigenval)-1),n_eigenval,1)
+ if(n_eigenval>=2){
+ for (i in 2:n_eigenval){
+ diagA$vectors<-cbind(diagA$vectors,
+ matrix(diagA$values[i]^(c(1:n_eigenval)-1),n_eigenval,1))
+ }
+ }
+ return(diagA)
+}
+
+
+StateVarX<-function(y,tt,X_0,B){
+ # The code obtains the first q-1 state variable using eq 5.1 in Brockwell, Davis and Yang 2011
+ Time<-length(tt)
+ q<-length(X_0)
+ e_q<-rep(0,q)
+ e_q[q]<-1
+ X_0<-matrix(X_0,q,1)
+ e_q<-matrix(e_q,q,1)
+ X<-matrix(0,q,Time)
+ int<-matrix(0,q,1)
+ for (i in c(2:Time)){
+ int<-int+(expm(B*(tt[i]-tt[(i-1)]))%*%(e_q*y[i]))*(tt[i]-tt[(i-1)])
+ X[i]<-as.matrix(expm(B*tt[i])%*%X_0+int)
+ }
+ return(X)
+}
+
+
+StateVarXp<-function(y,X_q,tt,B,q,p){
+ # The code computes the state variable X using the derivatives of eq 5.2
+ # see Brockwell, Davis and Yang 2011
+
+ diagMatB<-yuima.carma.eigen(B)
+ if(length(diagMatB$values)>1){
+ MatrD<-diag(diagMatB$values)
+ }else{
+ MatrD<-as.matrix(diagMatB$values)
+ }
+ MatrR<-diagMatB$vectors
+ idx.r<-c(q:(p-1))
+ elem.X <-length(idx.r)
+ YMatr<-matrix(0,q,length(y))
+ YMatr[q,]<-y
+ OutherX<-matrix(NA,q,length(y))
+ OutherXalt<-matrix(NA,q,length(y))
+ for(i in 1:elem.X){
+ OutherX[i,]<-((MatrR%*%MatrD^(idx.r[i])%*%solve(MatrR))%*%X_q+
+ (MatrR%*%MatrD^(idx.r[i]-1)%*%solve(MatrR))%*%YMatr)[1,]
+ }
+ X.StatVar<-rbind(X_q,OutherX)
+ return(X.StatVar)
+}
+
+bEvalPoly<-function(b,lambdax){
+ result<-sum(b*lambdax^(c(1:length(b))-1))
+ return(result)
+}
+
+aEvalPoly<-function(a,lambdax){
+ p<-length(a)
+ a.new<-c(1,a[c(1:(p-1))])
+ pa.new<-c(p:1)*a.new
+ result<-sum(pa.new*lambdax^(c(p:1)-1))
+ return(result)
+}
+
+
+CarmaRecovNoise<-function(yuima, param, data=NULL){
+ if( missing(param) )
+ yuima.stop("Parameter values are missing.")
+
+ if(!is.list(param))
+ yuima.stop("Argument 'param' must be of list type.")
+
+ vect.param<-as.numeric(param)
+ name.param<-names(param)
+ names(vect.param)<-name.param
+
+ if(is(yuima,"yuima")){
+ model<-yuima at model
+ if(is.null(data)){
+ observ<-yuima at data
+ }else{observ<-data}
+}else{
+ if(is(yuima,"yuima.carma")){
+ model<-yuima
+ if(is.null(data)){
+ yuima.stop("Missing data")
+ }
+ observ<-data
+ }
+}
+
+ if(!is(observ,"yuima.data")){
+ yuima.stop("Data must be an object of class yuima.data-class")
+ }
+
+ info<-model at info
+
+ numb.ar<-info at p
+ name.ar<-paste(info at ar.par,c(numb.ar:1),sep="")
+ ar.par<-vect.param[name.ar]
+
+ numb.ma<-info at q
+ name.ma<-paste(info at ma.par,c(0:numb.ma),sep="")
+ ma.par<-vect.param[name.ma]
+
+ loc.par=NULL
+ if (length(info at loc.par)!=0){
+ loc.par<-vect.param[info at loc.par]
+ }
+
+ scale.par=NULL
+ if (length(info at scale.par)!=0){
+ scale.par<-vect.param[info at scale.par]
+ }
+
+ lin.par=NULL
+ if (length(info at lin.par)!=0){
+ lin.par<-vect.param[info at lin.par]
+ }
+
+
+
+ ttt<-observ at zoo.data[[1]]
+ tt<-index(ttt)
+ y<-coredata(ttt)
+
+ levy<-yuima.CarmaRecovNoise(y,tt,ar.par,ma.par, loc.par, scale.par, lin.par)
+ inc.levy<-diff(t(levy))
+ return(inc.levy)
+}
+
+
+
+yuima.CarmaRecovNoise<-function(y,tt,ar.par,ma.par, loc.par=NULL, scale.par=NULL, lin.par=NULL){
+ if(!is.null(loc.par)){
+ yuima.warn("the loc.par will be implemented as soon as possible")
+ return(NULL)
+ }
+ if(!is.null(scale.par)){
+ yuima.warn("the scale.par will be implemented as soon as possible")
+ return(NULL)
+ }
+ if(!is.null(lin.par)){
+ yuima.warn("the lin.par will be implemented as soon as possible")
+ return(NULL)
+ }
+
+ p<-length(ar.par)
+ q<-length(ma.par)
+
+ if(q==1){
+ yuima.warn("the car(p) process will be implemented as soon as possible")
+ return(NULL)
+ }else{
+ A<-MatrixA(ar.par[c(p:1)])
+
+ b_q<-tail(ma.par,n=1)
+
+ newma.par<-ma.par/b_q
+ # We build the matrix B which is necessary for building eq 5.2
+ ynew<-y/b_q
+ B<-MatrixA(newma.par[c(1:(q-1))])
+ diagB<-yuima.carma.eigen(B)
+
+ e_q<-rep(0,(q-1))
+ if((q-1)>0){
+ e_q[(q-1)]<-1
+ X_0<-rep(0,(q-1))
+ }else{
+ e_q<-1
+ X_0<-0
+ }
+
+
+
+ X_q<-StateVarX(ynew,tt,X_0,B)
+
+ #plot(t(X_q))
+
+ X.StVa<-StateVarXp(ynew,X_q,tt,B,q-1,p)
+
+ #plot(y)
+
+ diagA<-yuima.carma.eigen(A)
+
+
+
+
+ BinLambda<-rep(NA,length(ma.par))
+ for(i in c(1:length(diagA$values))){
+ BinLambda[i]<-bEvalPoly(ma.par,diagA$values[i])
+ }
+ MatrBLam<-diag(BinLambda)
+
+ # We get the Canonical Vector Space in eq 2.17
+ Y_CVS<-MatrBLam%*%solve(diagA$vectors)%*%X.StVa #Canonical Vector Space
+
+ # We verify the prop 2 in the paper "Estimation for Non-Negative Levy Driven CARMA process
+ # yver<-Y_CVS[1,]+Y_CVS[2,]
+
+ # plot(yver)
+
+ # plot(y)
+
+
+ idx.r<-match(0,Im(diagA$values))
+ lambda.r<-diagA$values[idx.r]
+ int<-0
+ lev.und<-matrix(0,1,length(tt))
+ derA<-aEvalPoly(ar.par[c(p:1)],lambda.r)
+ for(t in c(2:length(tt))){
+ int<-int+y[t]*(tt[t]-tt[t-1])
+ lev.und[,t]<-derA/BinLambda[idx.r]*(Y_CVS[idx.r,t]-Y_CVS[idx.r,1]-lambda.r*int)
+ }
+ }
+ return(lev.und)
+}
+
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2014-01-12 16:34:58 UTC (rev 266)
+++ pkg/yuima/R/qmle.R 2014-01-13 21:07:48 UTC (rev 267)
@@ -112,23 +112,62 @@
}
if(is(yuima at model, "yuima.carma") && length(yuima at model@parameter at jump)!=0){
- 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)
+ CPlist <- c("dgamma", "dexp")
+ codelist <- c("rIG", "rgamma")
+ if(yuima at model@measure.type=="CP"){
+ 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,CPlist))){
+ yuima.warn("carma(p,q): the qmle for a carma(p,q) driven by a Compound Poisson with no-negative random size will be implemented as soon as possible")
+ return(NULL)
+ }
+ 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)
+ if(!is.na(match(measurefunc,codelist))){
+ yuima.warn("carma(p,q): the qmle for a carma(p,q) driven by a Compound Poisson with no-negative random size will be implemented as soon as possible")
+ 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(yuima at model, "yuima.carma") && 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
- jump.par <- yuima at model@parameter at jump
+#01/01 we introduce the new variable in order
+# to take into account the parameters in the starting conditions
+
+ if(is(yuima at model, "yuima.carma")){
+ #if(length(yuima at model@info at scale.par)!=0){
+ xinit.par <- yuima at model@parameter at xinit
+ #}
+ }
+
+
+ jump.par <- yuima at model@parameter at jump
measure.par <- yuima at model@parameter at measure
common.par <- yuima at model@parameter at common
-
- JointOptim <- joint
+
+ JointOptim <- joint
+ if(is(yuima at model, "yuima.carma") && length(yuima at model@parameter at jump)!=0){
+ if(any((match(jump.par, drift.par)))){
+ 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
@@ -137,7 +176,7 @@
# if(is(yuima at model, "yuima.carma")){
# JointOptim <- TRUE
# yuima.warm("Carma(p.q): The case of common parameters in Drift and Diffusion Term will be implemented as soon as possible,")
-# #return(null)
+# #return(NULL)
# }
}
@@ -158,6 +197,14 @@
if(length(drift.par)>0)
fullcoef <- c(fullcoef, drift.par)
+ if(is(yuima at model,"yuima.carma") &&
+ (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)
+
+ }
+
npar <- length(fullcoef)
fixed.par <- names(fixed)
@@ -167,6 +214,10 @@
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)]
@@ -174,7 +225,14 @@
idx.diff <- match(diff.par, nm)
idx.drift <- match(drift.par, nm)
- idx.fixed <- match(fixed.par, nm)
+ #01/01
+ if(is(yuima at model, "yuima.carma")){
+ # if(length(yuima at model@info at scale.par)!=0){
+ idx.xinit <- as.integer(na.omit(match(xinit.par,nm)))
+ #}
+ }
+
+ idx.fixed <- match(fixed.par, nm)
tmplower <- as.list( rep( -Inf, length(nm)))
names(tmplower) <- nm
@@ -335,11 +393,21 @@
mydots$upper <- unlist( upper[ nm[idx.drift] ])
mydots$lower <- unlist( lower[ nm[idx.drift] ])
if(length(mydots$par)>1){
- if(is(yuima at model, "yuima.carma")&& length(yuima at model@info at scale.par)!=0){
- name_b0<-paste0(yuima at model@info at ma.par,"0",collapse="")
- index_b0<-match(name_b0,nm)
- mydots$lower[index_b0]<-1
- mydots$upper[index_b0]<-1+10^(-7)
+ if(is(yuima at model, "yuima.carma")){
+ if(length(yuima at model@info at scale.par)!=0){
+ name_b0<-paste0(yuima at model@info at ma.par,"0",collapse="")
+ index_b0<-match(name_b0,nm)
+ mydots$lower[index_b0]<-1
+ mydots$upper[index_b0]<-1+10^(-7)
+ }
+ if (length(yuima at model@info at loc.par)!=0){
+ mydots$upper <- unlist( upper[ nm ])
+ mydots$lower <- unlist( lower[ nm ])
+ idx.tot<-unique(c(idx.drift,idx.xinit))
+ new.start <- start[idx.tot]
+ names(new.start) <- nm[idx.tot]
+ mydots$par <- unlist(new.start)
+ }
}
oout1 <- do.call(optim, args=mydots)
# oout1 <- optim(mydots$par,f,method = "L-BFGS-B" , lower = mydots$lower, upper = mydots$upper)
@@ -357,8 +425,14 @@
}
theta2 <- oout1$par
} ## endif(length(idx.drift)>0)
+
+
oout1 <- list(par= c(theta1, theta2))
- names(oout1$par) <- c(diff.par,drift.par)
+ if (! is(yuima at model,"yuima.carma")){
+ names(oout1$par) <- c(diff.par,drift.par)
+ }
+
+
oout <- oout1
} ### endif JointOptim
@@ -369,15 +443,19 @@
fDrift <- function(p) {
mycoef <- as.list(p)
- names(mycoef) <- drift.par
- mycoef[diff.par] <- coef[diff.par]
+ if(! is(yuima at model,"yuima.carma")){
+ names(mycoef) <- drift.par
+ mycoef[diff.par] <- coef[diff.par]
+ }
minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
}
fDiff <- function(p) {
- mycoef <- as.list(p)
- names(mycoef) <- diff.par
- mycoef[drift.par] <- coef[drift.par]
+ mycoef <- as.list(p)
+ if(! is(yuima at model,"yuima.carma")){
+ names(mycoef) <- diff.par
+ mycoef[drift.par] <- coef[drift.par]
+ }
minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
}
@@ -402,7 +480,21 @@
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)
-
+ if(is(yuima at model,"yuima.carma") && 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,
+ 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,
+ factr = 1e+07, pgtol = 0, tmax = 10, temp = 10)
+ }
# nmsC <- names(con)
# if (method == "Nelder-Mead")
# con$maxit <- 500
@@ -427,8 +519,13 @@
#
if(!HaveDriftHess & (length(drift.par)>0)){
#hess2 <- .Internal(optimhess(coef[drift.par], fDrift, NULL, conDrift))
- hess2 <- optimHess(coef[drift.par], fDrift, NULL, control=conDrift)
- HESS[drift.par,drift.par] <- hess2
+ if(!is(yuima at model,"yuima.carma")){
+ hess2 <- optimHess(coef[drift.par], fDrift, NULL, control=conDrift)
+ HESS[drift.par,drift.par] <- hess2
+ } else{
+ hess2 <- optimHess(coef, fDrift, NULL, control=conDrift)
+ HESS <- hess2
+ }
if(is(yuima at model,"yuima.carma") && 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))
@@ -461,9 +558,62 @@
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)
+# new("mle", call = call, coef = coef, fullcoef = unlist(mycoef),
+# vcov = vcov, min = min, details = oout, minuslogl = minusquasilogl,
+# method = method)
+#LM 11/01
+ final_res<-new("mle", call = call, coef = coef, fullcoef = unlist(mycoef),
+ vcov = vcov, min = min, details = oout, minuslogl = minusquasilogl,
+ method = method)
+
+if(!is(yuima at model,"yuima.carma")){
+ 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]
+ }
+
+ ttt<-observ at zoo.data[[1]]
+ tt<-index(ttt)
+ y<-coredata(ttt)
+
+ levy<-yuima.CarmaRecovNoise(y,tt,ar.par,ma.par, loc.par, scale.par, lin.par)
+ inc.levy<-NULL
+ if (!is.null(levy)){
+ inc.levy<-diff(t(levy))
+ }
+ carma_final_res<-list(mle=final_res,Incr=inc.levy,model=yuima)
+
+# the best should be to build a new class which extends both the mle and
+# the yuima class with an added slot that contains the estimated Levy increments
+
+
+ }
}
@@ -497,12 +647,22 @@
# 24/12
if(is(yuima at model, "yuima.carma") && 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
+# if(is(yuima at model, "yuima.carma")){
+# if(length(yuima at model@info at scale.par)!=0){
+# xinit.par <- yuima at model@parameter at xinit
+# }
+# }
+
+ if(is(yuima at model, "yuima.carma")){
+ xinit.par <- yuima at model@parameter at xinit
+ }
+
fullcoef <- NULL
if(length(diff.par)>0)
@@ -510,7 +670,19 @@
if(length(drift.par)>0)
fullcoef <- c(fullcoef, drift.par)
-
+ #01/01
+# if(is(yuima at model, "yuima.carma")){
+# if(length(yuima at model@info at scale.par)!=0){
+# if(length(xinit.par)>0)
+# fullcoef<-c(fullcoef, xinit.par)
+# }
+# }
+
+ if(is(yuima at model, "yuima.carma")){
+ if(length(xinit.par)>0)
+ fullcoef<-c(fullcoef, xinit.par)
+ }
+
# fullcoef <- c(diff.par, drift.par)
npar <- length(fullcoef)
# cat("\nparam\n")
@@ -526,18 +698,44 @@
idx.diff <- match(diff.par, nm)
idx.drift <- match(drift.par, nm)
-
+# if(is(yuima at model, "yuima.carma")){
+# if(length(yuima at model@info at scale.par)!=0){
+# idx.xinit <-as.integer(na.omit(match(xinit.par, nm)))
+# }
+# }
+
+ if(is(yuima at model, "yuima.carma")){
+ idx.xinit <-as.integer(na.omit(match(xinit.par, nm)))
+ }
+
h <- env$h
theta1 <- unlist(param[idx.diff])
theta2 <- unlist(param[idx.drift])
+
+
n.theta1 <- length(theta1)
n.theta2 <- length(theta2)
n.theta <- n.theta1+n.theta2
+#01/01
+# if(is(yuima at model, "yuima.carma")){
+# if(length(yuima at model@info at scale.par)!=0){
+# theta3 <- unlist(param[idx.xinit])
+# n.theta3 <- length(theta3)
+# n.theta <- n.theta1+n.theta2+n.theta3
+# }
+# }
- d.size <- yuima at model@equation.number
+ if(is(yuima at model, "yuima.carma")){
+ theta3 <- unlist(param[idx.xinit])
+ n.theta3 <- length(theta3)
+ n.theta <- n.theta1+n.theta2+n.theta3
+ }
+ d.size <- yuima at model@equation.number
+
+
n <- length(yuima)[1]
@@ -545,7 +743,7 @@
# 24/12
d.size <-1
# We build the two step procedure as described in
- if(length(yuima at model@info at scale.par)!=0){
+ # if(length(yuima at model@info at scale.par)!=0){
prova<-as.numeric(param)
names(prova)<-fullcoef[oo]
param<-prova[c(length(prova):1)]
@@ -553,25 +751,67 @@
y<-as.numeric(env$X)
tt<-env$time
p<-yuima at model@info at p
- a<-param[c(1:p)]
-# a_names<-names(param[c(1:p)])
-# names(a)<-a_names
- b<-param[c((p+1):(length(param)-p+1))]
-# b_names<-names(param[c((p+1):(length(param)-p+1))])
-# names(b)<-b_names
- if(length(b)==1){
- b<-1
- } else{
- indx_b0<-paste0(yuima at model@info at ma.par,"0",collapse="")
- b[indx_b0]<-1
- }
+ ar.par <- yuima at model@info at ar.par
+ name.ar<-paste0(ar.par, c(1:p))
+ q <- yuima at model@info at q
+ ma.par <- yuima at model@info at ma.par
+ name.ma<-paste0(ma.par, c(0:q))
+ if (length(yuima at model@info at loc.par)==0){
+
+ a<-param[name.ar]
+ # a_names<-names(param[c(1:p)])
+ # names(a)<-a_names
+ b<-param[name.ma]
+ # b_names<-names(param[c((p+1):(length(param)-p+1))])
+ # names(b)<-b_names
+ if(length(yuima at model@info at scale.par)!=0){
+ if(length(b)==1){
+ b<-1
+ } else{
+ indx_b0<-paste0(yuima at model@info at ma.par,"0",collapse="")
+ b[indx_b0]<-1
+ }
sigma<-tail(param,1)
+ }else {sigma<-1}
strLog<-yuima.carma.loglik1(y, tt, a, b, sigma)
+ }else{
+ # 01/01
+# ar.par <- yuima at model@info at ar.par
+# name.ar<-paste0(ar.par, c(1:p))
+ a<-param[name.ar]
+# ma.par <- yuima at model@info at ma.par
+# q <- yuima at model@info at q
+ name.ma<-paste0(ma.par, c(0:q))
+ b<-param[name.ma]
+ if(length(yuima at model@info at scale.par)!=0){
+ if(length(b)==1){
+ b<-1
+ } else{
+ indx_b0<-paste0(yuima at model@info at ma.par,"0",collapse="")
+ b[indx_b0]<-1
+ }
+ scale.par <- yuima at model@info at scale.par
+ sigma <- param[scale.par]
+ } else{sigma <- 1}
+ loc.par <- yuima at model@info at loc.par
+ mu <- param[loc.par]
+
+# if (length(b)==p){
+# # Be useful for carma driven by levy process
+# mean.y<-mu*tail(b,n=1)/tail(a,n=1)*sigma
+# }else{
+# mean.y<-0
+# }
+ y.start <- y-mu
+
+ strLog<-yuima.carma.loglik1(y.start, tt, a, b, sigma)
+ }
+
QL<-strLog$loglikCdiag
- }else {
- yuima.warn("carma(p,q): the scale parameter is equal to 1. We will implemented as soon as possible")
- return(NULL)
- }
+# }else {
+# yuima.warn("carma(p,q): the scale parameter is equal to 1. We will implemented as soon as possible")
+# return(NULL)
+# }
} else{
drift <- drift.term(yuima, param, env)
diff <- diffusion.term(yuima, param, env)
@@ -1088,15 +1328,19 @@
fDrift <- function(p) {
mycoef <- as.list(p)
- names(mycoef) <- drift.par
- mycoef[diff.par] <- coef[diff.par]
+ if(! is(yuima at model,"yuima.carma")){
+ names(mycoef) <- drift.par
+ mycoef[diff.par] <- coef[diff.par]
+ }
minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
}
fDiff <- function(p) {
mycoef <- as.list(p)
- names(mycoef) <- diff.par
- mycoef[drift.par] <- coef[drift.par]
+ if(! is(yuima at model,"yuima.carma")){
+ names(mycoef) <- diff.par
+ mycoef[drift.par] <- coef[drift.par]
+ }
minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
}
Modified: pkg/yuima/man/qmle.Rd
===================================================================
--- pkg/yuima/man/qmle.Rd 2014-01-12 16:34:58 UTC (rev 266)
+++ pkg/yuima/man/qmle.Rd 2014-01-13 21:07:48 UTC (rev 267)
@@ -46,7 +46,7 @@
% parameters. A function ml.pl estimates parameters of the sdeModel by
% maximizing the quasi-likelihood.
\code{qmle} behaves more likely the standard \code{mle} function in \pkg{stats4} and
- argument \code{method} is one of the methods available in \code{\link{optim}}.
+ argument \code{method} is one of the methods available in \code{\link{optim}}.
\code{lse} calculates least squares estimators of the drift parameters. This is
useful for initial guess of \code{qmle} estimation.
@@ -58,6 +58,12 @@
\value{
\item{QL}{a real value.}
\item{opt}{a list with components the same as 'optim' function.}
+ \item{carmaopt}{if the model is an object of \code{\link{yuima.carma-class}}, \code{qmle} returns a list containing:
+ \describe{
+ \item{mle}{contains an object of class \code{mle-class}.}
+ \item{Incr}{contains a vector of estimated noise.}
+ \item{model}{contains an object of class \code{\link{yuima-class}}.}
+ }}
}
\author{The YUIMA Project Team}
\note{
@@ -217,7 +223,7 @@
sigma=0.23))
)
-summary(carmaopt0)
+summary(carmaopt0$mle)
@@ -225,14 +231,12 @@
# carma(p=2,q=1) driven by a brownian motion without location parameter
mod1<-setCarma(p=2,
- q=1,
- scale.par="sigma")
+ q=1)
true.parm1 <-list(a1=1.39631,
a2=0.05029,
b0=1,
- b1=2,
- sigma=0.23)
+ b1=2)
samp1<-setSampling(Terminal=100,n=250)
set.seed(123)
@@ -242,13 +246,41 @@
system.time(
carmaopt1 <- qmle(sim1, start=list(a1=1.39631,a2=0.05029,
- b0=1,b1=2,
- sigma=0.23),joint=TRUE)
+ b0=1,b1=2),joint=TRUE)
)
-summary(carmaopt1)
+summary(carmaopt1$mle)
+
+# carma(p=2,q=1) driven by a compound poisson process where the jump size is normally distributed.
+
+mod2<-setCarma(p=2,
+ q=1,
+ measure=list(intensity="1",df=list("dnorm(z, 0, 1)")),
+ measure.type="CP")
+
+true.parm2 <-list(a1=1.39631,
+ a2=0.05029,
+ b0=1,
+ b1=2)
+
+samp2<-setSampling(Terminal=100,n=250)
+set.seed(123)
+sim2<-simulate(mod2,
+ true.parameter=true.parm2,
+ sampling=samp2)
+
+system.time(
+ carmaopt2 <- qmle(sim2, start=list(a1=1.39631,a2=0.05029,
+ b0=1,b1=2),joint=TRUE)
+)
+
+summary(carmaopt2$mle)
+
+
+
}
+
% Add one or more standard keywords, see file 'KEYWORDS' in the
% R documentation directory.
\keyword{ts}
More information about the Yuima-commits
mailing list