[Yuima-commits] r721 - in pkg/yuima: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Feb 19 20:09:03 CET 2020
Author: lorenzo
Date: 2020-02-19 20:09:03 +0100 (Wed, 19 Feb 2020)
New Revision: 721
Modified:
pkg/yuima/R/qmleLevy.R
pkg/yuima/man/qmleLevy.Rd
Log:
updated qmleLevy.R and qmleLevy.Rd
Modified: pkg/yuima/R/qmleLevy.R
===================================================================
--- pkg/yuima/R/qmleLevy.R 2020-02-06 21:53:03 UTC (rev 720)
+++ pkg/yuima/R/qmleLevy.R 2020-02-19 19:09:03 UTC (rev 721)
@@ -3,8 +3,15 @@
########################################################################
-qmleLevy<-function(yuima,start,lower,upper,joint = FALSE,third = FALSE)
+qmleLevy<-function(yuima,start,lower,upper,joint = FALSE,third = FALSE,
+ Est.Incr = c("NoIncr","Incr","IncrPar"),
+ aggregation = TRUE)
{
+ cat("\nStarting QGMLE for SDE ... \n")
+ parameter<-yuima at model@parameter at all
+ orig.mylaw<-yuima at model@measure
+ mylaw<-yuima at model@measure$df
+
if(missing(yuima))
yuima.stop("yuima object is missing.")
@@ -11,6 +18,20 @@
if(!is(yuima,"yuima"))
yuima.stop("This function is for yuima-class.")
+ if(length(yuima at model@parameter at jump)>0)
+ paracoef <- yuima at model@parameter at jump
+
+ if(length(yuima at model@parameter at drift)>0)
+ paracoef <- c(paracoef, yuima at model@parameter at drift)
+ if(Est.Incr == "IncrPar"){
+ start0<-start
+ lower0<-lower
+ upper0<-upper
+ lev.names<-yuima at model@parameter at measure
+ }
+
+ DRIFT <- yuima at model@drift
+ JUMP <- yuima at model@jump.coeff
sdeModel<-yuima at model
if(length(sdeModel at parameter@measure)!=0){
nPar<-length(sdeModel at parameter@measure)
@@ -269,10 +290,137 @@
yuima at model@measure.type <- character(0)
yuima at model@parameter at jump <- character(0)
yuima at model@parameter at measure <- character(0)
- jres<-qmle(yuima,start = start,lower = lower,upper = upper,rcpp = TRUE, joint = TRUE,method = "L-BFGS-B")
- res<-list(joint = jres at coef)
+ # jres<-qmle(yuima,start = start,lower = lower,upper = upper,rcpp = TRUE, joint = TRUE,method = "L-BFGS-B")
+ # res<-list(joint = jres at coef)
+ res<-qmle(yuima,start = start,lower = lower,upper = upper,rcpp = TRUE, joint = TRUE,
+ method = "L-BFGS-B")
}
- res
+ if(Est.Incr == "NoIncr"){
+ return(res)
+ }
+ cat("\nStarting Estimation of Levy Increments ... \n")
+ data <- get.zoo.data(yuima)
+ s.size<-yuima at sampling@n
+ X<-as.numeric(data[[1]])
+ pX<-X[1:(s.size-1)]
+ inc<-double(s.size-1)
+ inc<-X[2:(s.size)]-pX
+ modeltime<-yuima at model@time.variable
+ modelstate<-yuima at model@solve.variable
+ tmp.env<-new.env()
+
+ # DRIFT <- yuima at model@drift
+ # JUMP <- yuima at model@jump.coeff
+
+ if(length(yuima at model@solve.variable)==1){
+ #parameter<-yuima at model@parameter at all
+ if(joint){
+ coeffic<- coef(res)
+ }else{
+ coeffic<- res[[1]]
+ if(length(res)>1){
+ for(j in c(2:length(res))){
+ coeffic<-c(coeffic,res[[j]])
+ }
+ }
+
+ }
+ mp<-match(names(coeffic),parameter)
+ esort <- coeffic[order(mp)]
+ for(i in 1:length(coeffic))
+ {
+ assign(parameter[i],esort[[i]],envir=tmp.env)
+ }
+
+ resi<-double(s.size-1)
+ assign(modeltime,yuima at sampling@delta,envir=tmp.env)
+ h<-yuima at sampling@delta
+ assign(modelstate,pX,envir=tmp.env)
+ jump.term<-eval(JUMP[[1]],envir=tmp.env)
+ drif.term<-eval(DRIFT,envir=tmp.env)
+ if(length(jump.term)==1){
+ jump.term <- rep(jump.term, s.size)
+ }
+ if(length(drif.term)==1){
+ drif.term <- rep(drif.term, s.size)
+ } # vectorization (note. if an expression type object does not include state.variable, the length of the item after "eval" operation is 1.)
+ for(s in 1:(s.size-1)){
+ nova<-sqrt((jump.term)^2) # normalized variance
+ resi[s]<-(1/(nova[s]))*(inc[s]-h*drif.term[s])
+ }
+ if(aggregation){
+ Ter <- yuima at sampling@Terminal
+ ures <- numeric(floor(Ter))
+ for(l in 1:floor(Ter)){
+ ures[l] <- sum(resi[(floor((l-1)/h)):(floor(l/h)-1)])
+ }
+ res.incr<-ures
+ }else{
+ res.incr<-resi
+ }
+
+ }else{}
+
+ if(Est.Incr == "Incr"){
+ return(list(res=res,Est.Incr=res.incr))
+ }
+ cat("\nEstimation Levy parameters ... \n")
+
+ if(class(mylaw)=="yuima.law"){
+ if(aggregation){
+ minusloglik <- function(para){
+ para[length(para)+1]<-1
+ names(para)[length(para)]<-yuima at model@time.variable
+ -sum(dens(mylaw, res.incr, param = para, log = T),
+ na.rm = T)
+ }
+ }else{
+ minusloglik <- function(para){
+ para[length(para)+1] <- yuima at sampling@delta
+ names(para)[length(para)]<-yuima at model@time.variable
+ -sum(dens(mylaw, res.incr, param = para, log = T),
+ na.rm = T)
+ }
+ }
+ para <- start0[lev.names]
+ lowerjump <- lower0[lev.names]
+ upperjump <- upper0[lev.names]
+ esti <- optim(fn = minusloglik, lower = lowerjump, upper = upperjump,
+ par = para, method = "L-BFGS-B")
+ return(list(res=res,Est.Incr=res.incr, meas=esti$par))
+ }else{
+ dist <- substr(as.character(orig.mylaw$df$expr), 2, 10^3)
+
+ startjump <- start0[lev.names]
+ lowerjump <- lower0[lev.names]
+ upperjump <- upper0[lev.names]
+
+ if(length(startjump) == 1){
+ logdens <- function(para){
+ exlogdens <- parse(text = sprintf("log(d%s)", dist))
+ assign(yuima at model@jump.variable, ures, envir = tmp.env)
+ assign(yuima at model@parameter at measure, para, envir = tmp.env)
+ sum(eval(exlogdens, envir = tmp.env))
+ }
+ intervaljump <- c(lowerjump[[1]], upperjump[[1]])
+ esti <- optimize(logdens, interval = intervaljump, maximum = TRUE)
+ return(list(sde=esort, meas=esti$maximum))
+ }else{
+ logdens <- function(para){
+ exlogdens <- parse(text = sprintf("log(d%s)", dist))
+ assign(yuima at model@jump.variable, ures, envir = tmp.env)
+ for(i in c(1:length(yuima at model@parameter at measure)))
+ assign(yuima at model@parameter at measure[i], para[[yuima at model@parameter at measure[i]]], envir = tmp.env)
+
+ sum(eval(exlogdens, envir = tmp.env))
+ }
+
+ esti <- optim(fn=logdens, lower = lowerjump, upper = upperjump, par = startjump,
+ method = "L-BFGS-B", control = list(fnscale = -1))
+ return(list(sde=esort, meas=esti$par))
+ }
+ }
+
}
Modified: pkg/yuima/man/qmleLevy.Rd
===================================================================
--- pkg/yuima/man/qmleLevy.Rd 2020-02-06 21:53:03 UTC (rev 720)
+++ pkg/yuima/man/qmleLevy.Rd 2020-02-19 19:09:03 UTC (rev 721)
@@ -1,6 +1,8 @@
\encoding{UTF-8}
\name{qmleLevy}
\alias{qmleLevy}
+\alias{Estimation.LevyIncr}
+\alias{LevySDE}
%- Also NEED an '\alias' for EACH other topic documented here.
\title{
Gaussian quasi-likelihood estimation for Levy driven SDE
@@ -9,7 +11,9 @@
Calculate the Gaussian quasi-likelihood and Gaussian quasi-likelihood estimators of Levy driven SDE.
}
\usage{
-qmleLevy(yuima, start, lower, upper, joint = FALSE, third = FALSE)
+qmleLevy(yuima, start, lower, upper, joint = FALSE,
+third = FALSE, Est.Incr = c("NoIncr", "Incr", "IncrPar"),
+aggregation = TRUE)
}
\arguments{
\item{yuima}{a yuima object.}
@@ -16,9 +20,11 @@
\item{lower}{a named list for specifying lower bounds of parameters.}
\item{upper}{a named list for specifying upper bounds of parameters.}
\item{start}{initial values to be passed to the optimizer.}
- \item{joint}{perform joint estimation or two stage estimation? by default joint=FALSE. If there exists an overlapping parameter, joint=TRUE does not work for the theoretical reason}
- \item{third}{perform third estimation? by default third=FALSE. If there exists an overlapping parameter, third=TRUE does not work for the
+ \item{joint}{perform joint estimation or two stage estimation, by default \code{joint=FALSE}. If there exists an overlapping parameter, \code{joint=TRUE} does not work for the theoretical reason}
+ \item{third}{perform third estimation by default \code{third=FALSE}. If there exists an overlapping parameter, \code{third=TRUE} does not work for the
theoretical reason.}
+ \item{Est.Incr}{the qmleLevy returns an object of \code{mle-clas}, by default \code{Est.Incr="NoIncr"}.}
+ \item{aggregation}{If \code{aggregation=TRUE}, the function returns the unit-time Levy increments. If \code{Est.Incr="IncrPar"}, the function estimates Levy parameters using the unit-time Levy increments.}
}
\details{
This function performs Gaussian quasi-likelihood estimation for Levy driven SDE.
More information about the Yuima-commits
mailing list