[Yuima-commits] r270 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jan 27 20:59:43 CET 2014
Author: lorenzo
Date: 2014-01-27 20:59:43 +0100 (Mon, 27 Jan 2014)
New Revision: 270
Added:
pkg/yuima/man/yuima.carma.qmle-class.Rd
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/R/AllClasses.R
pkg/yuima/R/qmle.R
pkg/yuima/R/sim.euler.R
pkg/yuima/man/CarmaRecovNoise.Rd
pkg/yuima/man/qmle.Rd
pkg/yuima/man/yuima.carma-class.Rd
Log:
Added new class yuima.carma.qmle and documentation
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2014-01-22 12:56:16 UTC (rev 269)
+++ pkg/yuima/DESCRIPTION 2014-01-27 19:59:43 UTC (rev 270)
@@ -3,7 +3,7 @@
Title: The YUIMA Project package (unstable version)
Version: 0.1.221
Date: 2014-01-13
-Depends: methods, zoo, stats4, utils, Matrix
+Depends: methods, zoo, stats4, utils, Matrix, GeneralizedHyperbolic
Suggests: cubature, mvtnorm
Author: YUIMA Project Team.
Maintainer: Stefano M. Iacus <stefano.iacus at R-project.org>
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2014-01-22 12:56:16 UTC (rev 269)
+++ pkg/yuima/NAMESPACE 2014-01-27 19:59:43 UTC (rev 270)
@@ -6,6 +6,7 @@
importFrom("Matrix")
importFrom(stats, confint)
importFrom(stats4)
+importFrom("GeneralizedHyperbolic")
importFrom(utils, toLatex)
@@ -22,7 +23,8 @@
"yuima.model",
"model.parameter",
"yuima.carma",
- "carma.info"
+ "carma.info",
+ "yuima.carma.qmle"
)
exportMethods(
Modified: pkg/yuima/R/AllClasses.R
===================================================================
--- pkg/yuima/R/AllClasses.R 2014-01-22 12:56:16 UTC (rev 269)
+++ pkg/yuima/R/AllClasses.R 2014-01-27 19:59:43 UTC (rev 270)
@@ -120,3 +120,12 @@
functional = "yuima.functional"
)
)
+
+# Class yuima.carma.qmle
+setClass("yuima.carma.qmle",representation(Incr.Lev = "matrix",
+ model = "yuima.carma"
+ ),
+ contains="mle"
+ )
+# The yuima.carma.qmle extends the S4 class "mle". It contains three slots: Estimated Levy,
+# The description of the carma model and the mle.
\ No newline at end of file
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2014-01-22 12:56:16 UTC (rev 269)
+++ pkg/yuima/R/qmle.R 2014-01-27 19:59:43 UTC (rev 270)
@@ -74,7 +74,7 @@
qmle <- function(yuima, start, method="BFGS", fixed = list(), print=FALSE,
- lower, upper, joint=FALSE, ...){
+ lower, upper, joint=FALSE, Est.Incr=TRUE, ...){
if(is(yuima at model, "yuima.carma")){
NoNeg.Noise<-FALSE
}
@@ -128,20 +128,21 @@
}
# 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 non-Negative Levy will be implemented as soon as possible")
- NoNeg.Noise<-TRUE
- if((yuima at model@info at q+1)==(yuima at model@info at q+1)){
- start[[mean.noise]]<-1
- }
- #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 non-Negative Levy will be implemented as soon as possible")
+ NoNeg.Noise<-TRUE
+ if((yuima at model@info at q+1)==(yuima at model@info at q+1)){
+ 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 ")
@@ -624,9 +625,21 @@
# 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")){
+ final_res<-new("mle", call = call, coef = coef, fullcoef = unlist(mycoef),
+ vcov = vcov, min = min, details = oout, minuslogl = minusquasilogl,
+ method = method)
+ }else{
+ if(Est.Incr==TRUE){
+ final_res<-new("yuima.carma.qmle", call = call, coef = coef, fullcoef = unlist(mycoef),
+ vcov = vcov, min = min, details = oout, minuslogl = minusquasilogl,
+ method = method)
+ }else{
+ 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)
@@ -670,12 +683,155 @@
if (!is.null(levy)){
inc.levy<-diff(t(levy))
}
- carma_final_res<-list(mle=final_res,Incr=inc.levy,model=yuima)
+ # INSERT HERE THE NECESSARY STEPS FOR FINDING THE PARAMETERS OF LEVY
-# 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
+ 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)
+ names.par<-c(unique(c(drift.par,diff.par)),unique(c(measure.par)))
+ 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))
+ rownames(cov)<-names.par
+ colnames(cov)<-names.par
+ if(is.null(loc.par)){
+ cov[unique(c(drift.par,diff.par)),unique(c(drift.par,diff.par))]<-dummycovCarmapar
+ }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"){
+ # tmp <- regexpr("\\(", yuima at model@measure$df$exp)[1]
+ # measurefunc <- substring(yuima at model@measure$df$exp, 1, tmp-1)
+ if(measurefunc=="dnorm"){
+
+ }
+ if(measurefunc=="dgamma"){
+
+ }
+ if(measurefunc=="dexp"){
+
+ }
+ }
+ if(yuima at model@measure.type=="code"){
+ # # "rIG", "rNIG", "rgamma", "rbgamma", "rngamma"
+
+
+ if(measurefunc=="rIG"){
+ # result.Levy<-gigFit(inc.levy)
+ # Inc.Parm<-coef(result.Levy)
+ # IncVCOV<--solve(gigHessian(inc.levy, param=Inc.Parm))
+ }
+ if(measurefunc=="rNIG"){
+ # inc.levy
+ # measure.par
+ # sapply(gregexpr("\\W+", measurefunc),length)
+
+ name.func.dummy <- as.character(model at measure$df$expr[1])
+ name.func<- substr(name.func.dummy,1,(nchar(name.func.dummy)-1))
+ names.measpar<-rev(as.vector(strsplit(name.func,', '))[[1]][-1])
+ valuemeasure<-as.numeric(names.measpar)
+ NaIdx<-which(is.na(valuemeasure))
+ if(length(NaIdx)!=0){
+ yuima.warn("the constrained MLE for levy increment will be implemented as soon as possible")
+ }
+
+ inc.levy1<-diff(cumsum(inc.levy)[seq(from=1,
+ to=yuima at sampling@n[1],
+ by=(yuima at sampling@n/yuima at sampling@Terminal)[1]
+ )])
+ result.Levy<-nigFit(inc.levy1)
+
+ Inc.Parm<-coef(result.Levy)
+ IncVCOV<--solve(nigHessian(inc.levy, param=Inc.Parm))
+
+ names(Inc.Parm)[NaIdx]<-measure.par
+ #prova<-as.matrix(IncVCOV)
+ rownames(IncVCOV)[NaIdx]<-as.character(measure.par)
+ colnames(IncVCOV)[NaIdx]<-as.character(measure.par)
+
+
+ coef<-NULL
+ coef<-c(dummycoeffCarmapar,Inc.Parm)
+ # names.par<-c(unique(c(drift.par,diff.par)),names(Inc.Parm))
+ #
+ names.par<-names(coef)
+ cov<-NULL
+ cov<-matrix(0,length(names.par),length(names.par))
+ rownames(cov)<-names.par
+ colnames(cov)<-names.par
+ if(is.null(loc.par)){
+ cov[unique(c(drift.par,diff.par)),unique(c(drift.par,diff.par))]<-dummycovCarmapar
+ }else{
+ 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(measurefunc=="rgamma"){
+ # result.Levy<-gigFit(inc.levy)
+ # Inc.Parm<-coef(result.Levy)
+ # IncVCOV<--solve(gigHessian(inc.levy, param=Inc.Parm))
+
+ }
+ if(measurefunc=="rbgamma"){
+
+ }
+ if(measurefunc=="rngamma"){
+
+ }
+
+
+
+
+ }
+ }
+# dummycovCarmapar<-vcov[unique(c(drift.par,diff.par)),unique(c(drift.par,diff.par))]
+# dummycovCarmaNoise<-vcov[unique(measure.par),unique(c(measure.par))] #we need to adjusted
+# dummycoeffCarmapar<-coef[unique(c(drift.par,diff.par))]
+# dummycoeffCarmaNoise<-coef[unique(c(measure.par))]
+# coef<-NULL
+# coef<-c(dummycoeffCarmapar,dummycoeffCarmaNoise)
+# names.par<-c(unique(c(drift.par,diff.par)),unique(c(measure.par)))
+# names(coef)<-names.par
+# cov<-NULL
+# cov<-matrix(0,length(names.par),length(names.par))
+# rownames(cov)<-names.par
+# colnames(cov)<-names.par
+# cov[unique(c(drift.par,diff.par)),unique(c(drift.par,diff.par))]<-dummycovCarmapar
+# cov[unique(c(measure.par)),unique(c(measure.par))]<-dummycovCarmaNoise
+
+# carma_final_res<-list(mle=final_res,Incr=inc.levy,model=yuima)
+ if(Est.Incr==TRUE){
+ # START FROM HERE 24/01
+ carma_final_res<-new("yuima.carma.qmle", call = call, coef = coef, fullcoef = unlist(mycoef),
+ vcov = cov, min = min, details = oout, minuslogl = minusquasilogl,
+ method = method, Incr.Lev = inc.levy,
+ model = yuima at model)
+ }else{
+ carma_final_res<-new("mle", call = call, coef = coef, fullcoef = unlist(mycoef),
+ vcov = cov, min = min, details = oout, minuslogl = minusquasilogl,
+ method = method)
+ }
+ return(carma_final_res)
}
}
@@ -909,7 +1065,9 @@
if (length(b)==p){
mean.noise<-param[mean.noise]
# Be useful for carma driven by levy process
- mean.y<-mean.noise*tail(b,n=1)/tail(a,n=1)*sigma
+ # mean.y<-mean.noise*tail(b,n=1)/tail(a,n=1)*sigma
+ mean.y<-mean(y-mu)
+
}else{
mean.y<-0
}
@@ -1531,4 +1689,10 @@
method = method)
}
+# Plot Method for yuima.carma.qmle
+setMethod("plot",signature(x="yuima.carma.qmle"),
+ function(x, ...){
+ plot(x at Incr.Lev, ...)
+ }
+)
Modified: pkg/yuima/R/sim.euler.R
===================================================================
--- pkg/yuima/R/sim.euler.R 2014-01-22 12:56:16 UTC (rev 269)
+++ pkg/yuima/R/sim.euler.R 2014-01-27 19:59:43 UTC (rev 270)
@@ -163,6 +163,17 @@
}
#lambda <- integrate(sdeModel at measure$df$func, 0, Inf)$value * eta0
#lambda <- integrate(tmp.expr, 0, Inf)$value * eta0 ##bug:2013/10/28
+
+ dummyList<-as.list(env)
+ lgth.meas<-length(yuima at model@parameter at measure)
+ if(lgth.meas>1){
+ for(i in c(2:lgth.meas)){
+ idx.dummy<-yuima at model@parameter at measure[i]
+ assign(idx.dummy,as.numeric(dummyList[idx.dummy]))
+ }
+ }
+
+
lambda <- integrate(tmp.expr, -Inf, Inf)$value * eta0
##:: lambda = nu() (p6)
@@ -233,6 +244,14 @@
## rstable=paste("rstable(n, ", args[2], ", ", args[3], ", ", args[4], ", ", args[5], ", ", args[6], ")")
rstable=paste("rstable(n, ", args[2], ", ", args[3], ", ", args[4], "*delta^(1/",args[2],"), ", args[5], "*delta)")
)
+ dummyList<-as.list(env)
+ lgth.meas<-length(yuima at model@parameter at measure)
+ if(lgth.meas!=0){
+ for(i in c(1:lgth.meas)){
+ idx.dummy<-yuima at model@parameter at measure[i]
+ assign(idx.dummy,as.numeric(dummyList[idx.dummy]))
+ }
+ }
if(is.null(dZ)){ ##:: "otherwise"
cat(paste("Code \"", code, "\" not supported yet.\n", sep=""))
Modified: pkg/yuima/man/CarmaRecovNoise.Rd
===================================================================
--- pkg/yuima/man/CarmaRecovNoise.Rd 2014-01-22 12:56:16 UTC (rev 269)
+++ pkg/yuima/man/CarmaRecovNoise.Rd 2014-01-27 19:59:43 UTC (rev 270)
@@ -87,10 +87,9 @@
# We estimate the parameter using qmle.
carmaopt1 <- qmle(sim1, start=true.parm1)
-summary(carmaopt1$mle)
+summary(carmaopt1)
# Internally qmle uses CarmaRecovNoise. The result is in
-carmaopt1$Incr->inc.Levy1
-plot(inc.Levy1)
+plot(carmaopt1)
# Ex.3: Carma(p=2,q=1) with scale and location parameters
# driven by a Compound Poisson
@@ -118,11 +117,10 @@
# We estimate the Carma and we plot the underlying noise.
carmaopt2 <- qmle(sim2, start=true.parm2)
-summary(carmaopt2$mle)
+summary(carmaopt2)
-carmaopt2$Incr->inc.Levy2
# Increments estimated by CarmaRecovNoise
-plot(inc.Levy2)
+plot(carmaopt2)
}
Modified: pkg/yuima/man/qmle.Rd
===================================================================
--- pkg/yuima/man/qmle.Rd 2014-01-22 12:56:16 UTC (rev 269)
+++ pkg/yuima/man/qmle.Rd 2014-01-27 19:59:43 UTC (rev 270)
@@ -19,7 +19,7 @@
%ml.ql(yuima,theta2,theta1,h,theta2.lim=matrix(c(0,1),1,2),theta1.lim=matrix(c(0,1),1,2),print=FALSE,method,param,interval)
%ql(yuima,theta2,theta1,h,print=FALSE,param)
%rql(yuima,theta2,theta1,ptheta2,ptheta1,h,print=FALSE,param,prevparam)
-qmle(yuima, start, method="BFGS", fixed = list(), print=FALSE, lower, upper, joint=FALSE, ...)
+qmle(yuima, start, method="BFGS", fixed = list(), print=FALSE, lower, upper, joint=FALSE, Est.Incr=TRUE, ...)
quasilogl(yuima, param, print=FALSE)
}
\arguments{
@@ -39,6 +39,7 @@
\item{start}{initial values to be passed to the optimizer.}
\item{fixed}{for conditional (quasi)maximum likelihood estimation.}
\item{joint}{perform joint estimation or two stage estimation? by default \code{joint=FALSE}.}
+ \item{Est.Incr}{If the yuima model is an object of \code{\link{yuima.carma-class}} the \code{qmle} returns an object of \code{\link{yuima.carma.qmle-class}} or object of class \code{mle-class}. By default \code{Est.Incr=TRUE}.}
\item{...}{passed to \code{\link{optim}} method. See Examples.}
}
\details{
@@ -58,12 +59,7 @@
\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}}.}
- }}
+ \item{carmaopt}{if the model is an object of \code{\link{yuima.carma-class}}, \code{qmle} returns an object \code{\link{yuima.carma.qmle-class}}}
}
\author{The YUIMA Project Team}
\note{
@@ -225,11 +221,12 @@
sigma=0.23))
)
-summary(carmaopt0$mle)
+summary(carmaopt0)
+
# carma(p=2,q=1) driven by a brownian motion without location parameter
mod1<-setCarma(p=2,
@@ -251,8 +248,10 @@
b0=1,b1=2),joint=TRUE)
)
-summary(carmaopt1$mle)
+summary(carmaopt1)
+plot(carmaopt1)
+
# carma(p=2,q=1) driven by a compound poisson process where the jump size is normally distributed.
mod2<-setCarma(p=2,
@@ -276,10 +275,41 @@
b0=1,b1=2),joint=TRUE)
)
-summary(carmaopt2$mle)
+summary(carmaopt2)
+plot(carmaopt2)
+# carma(p=2,q=1) driven by a normal inverse gaussian process
+mod3<-setCarma(p=2,q=1,
+ measure=list(df=list("rNIG(z, alpha, beta, delta, mu)")),
+ measure.type="code")
+#
+# True param
+true.param3<-list(a1=1.39631,
+ a2=0.05029,
+ b0=1,
+ b1=2,
+ alpha=1,
+ beta=0,
+ delta=1,
+ mu=0)
+
+samp3<-setSampling(Terminal=100,n=200)
+set.seed(123)
+
+sim3<-simulate(mod3,
+ true.parameter=true.param3,
+ sampling=samp3)
+
+
+carmaopt3<-qmle(sim3,start=true.param3)
+
+summary(carmaopt3)
+
+plot(carmaopt3)
+
+
}
Modified: pkg/yuima/man/yuima.carma-class.Rd
===================================================================
--- pkg/yuima/man/yuima.carma-class.Rd 2014-01-22 12:56:16 UTC (rev 269)
+++ pkg/yuima/man/yuima.carma-class.Rd 2014-01-27 19:59:43 UTC (rev 270)
@@ -54,6 +54,7 @@
\code{\link{simulate}}.}
\item{toLatex}{This method converts an object of \code{yuima.carma-class} to character vectors with LaTeX markup.}
\item{CarmaRecovNoise}{Recovering underlying Levy. For more information see \code{\link{CarmaRecovNoise}}. }
+ \item{qmle}{Quasi maximum likelihood estimation procedure. For more information see \code{\link{qmle}}. }
}
}
\author{The YUIMA Project Team}
Added: pkg/yuima/man/yuima.carma.qmle-class.Rd
===================================================================
--- pkg/yuima/man/yuima.carma.qmle-class.Rd (rev 0)
+++ pkg/yuima/man/yuima.carma.qmle-class.Rd 2014-01-27 19:59:43 UTC (rev 270)
@@ -0,0 +1,33 @@
+\name{yuima.carma.qmle-class}
+\docType{class}
+\alias{yuima.carma.qmle-class}
+\alias{plot,yuima.carma.qmle,ANY-method}
+\alias{qmle.carma}
+\alias{carma.qmle}
+%%\alias{setSampling,yuima.carma-method}
+
+\title{Class for Quasi Maximum Likelihood Estimation of CARMA(p,q) model}
+\description{
+ The \code{yuima.carma.qmle} class is a class of the \pkg{yuima} package that extends the \code{mle-class} of the \pkg{stats4} package.
+}
+\section{Slots}{
+ \describe{
+ \item{\code{Incr.Lev}:}{is an \code{matrix} object that contains the estimated increments of the noise obtained using \code{\link{CarmaRecovNoise}}.}
+ \item{\code{model}:}{is an object of of \code{\link{yuima.carma-class}}.}
+ \item{\code{call}:}{is an object of class \code{language}. }
+ \item{\code{coef}:}{is an object of class \code{numeric} that contains estimated parameters.}
+ \item{\code{fullcoef}:}{is an object of class \code{numeric} that contains estimated and fixed parameters.}
+ \item{\code{vcov}:}{is an object of class \code{matrix}.}
+ \item{\code{min}:}{is an object of class \code{numeric}.}
+ \item{\code{minuslogl}:}{is an object of class \code{function}.}
+ \item{\code{method}:}{is an object of class \code{character}.}
+ }
+}
+\section{Methods}{
+ \describe{
+ \item{plot}{Plot method for estimated increment of the noise.}
+ \item{Methods mle}{All methods for \code{mle-class} are available.}
+ }
+}
+\author{The YUIMA Project Team}
+\keyword{classes}
More information about the Yuima-commits
mailing list