[Yuima-commits] r766 - in pkg/yuima: . R man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Dec 15 14:53:54 CET 2021
Author: phoenix844
Date: 2021-12-15 14:53:54 +0100 (Wed, 15 Dec 2021)
New Revision: 766
Added:
pkg/yuima/R/ae.R
pkg/yuima/man/ae.Rd
pkg/yuima/man/aeCharacteristic.Rd
pkg/yuima/man/aeDensity.Rd
pkg/yuima/man/aeExpectation.Rd
pkg/yuima/man/aeKurtosis.Rd
pkg/yuima/man/aeMarginal.Rd
pkg/yuima/man/aeMean.Rd
pkg/yuima/man/aeMoment.Rd
pkg/yuima/man/aeSd.Rd
pkg/yuima/man/aeSkewness.Rd
pkg/yuima/man/yuima.ae-class.Rd
pkg/yuima/src/ae.cpp
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/NEWS
pkg/yuima/src/yuima_init.c
Log:
Add asymptotic expansion
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2021-11-28 14:19:07 UTC (rev 765)
+++ pkg/yuima/DESCRIPTION 2021-12-15 13:53:54 UTC (rev 766)
@@ -1,11 +1,11 @@
Package: yuima
Type: Package
Title: The YUIMA Project Package for SDEs
-Version: 1.11.0
+Version: 1.12.0
Depends: R(>= 2.10.0), methods, zoo, stats4, utils, expm, cubature,
mvtnorm
Imports: Rcpp (>= 0.12.1), boot (>= 1.3-2), glassoFast, wavethresh,
- coda
+ coda, calculus (>= 0.2.0)
Author: YUIMA Project Team
Maintainer: Stefano M. Iacus <stefano.iacus at unimi.it>
Description: Simulation and Inference for SDEs and Other Stochastic Processes.
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2021-11-28 14:19:07 UTC (rev 765)
+++ pkg/yuima/NAMESPACE 2021-12-15 13:53:54 UTC (rev 766)
@@ -68,6 +68,16 @@
importFrom(coda, as.mcmc)
importFrom("stats", "cov", "var")
+# added by EG on Dec. 15, 2021
+importFrom(calculus,"%dot%")
+importFrom(calculus,"%gradient%")
+importFrom(calculus,"%inner%")
+importFrom(calculus,"%mx%")
+importFrom(calculus,"%outer%")
+importFrom(calculus,"%prod%")
+importFrom(calculus,"%sum%")
+
+
exportClasses("yuima",
"yuima.data",
"yuima.law",
@@ -96,7 +106,8 @@
"yuima.PPR",
"yuima.Hawkes",
"yuima.PPR.qmle",
-"yuima.snr"
+"yuima.snr",
+"yuima.ae" # added by EG on Dec. 15, 2021
)
exportMethods(
@@ -231,6 +242,18 @@
export(Diagnostic.Cogarch)
export(Diagnostic.Carma)
+# Asymptotic expansion (added by EG on Dec. 15, 2021)
+export(ae)
+export(aeCharacteristic)
+export(aeDensity)
+export(aeExpectation)
+export(aeKurtosis)
+export(aeMarginal)
+export(aeMean)
+export(aeMoment)
+export(aeSd)
+export(aeSkewness)
+
# Methods
export(rand)# random number generator of a Levy process specified by user
export(dens)
Modified: pkg/yuima/NEWS
===================================================================
--- pkg/yuima/NEWS 2021-11-28 14:19:07 UTC (rev 765)
+++ pkg/yuima/NEWS 2021-12-15 13:53:54 UTC (rev 766)
@@ -79,4 +79,7 @@
2021/9/24: fixed a bug regarding the "fixed" argument in qmle
2021/11/23: fixed a bug regarding the "psd" argument in lmm
2021/11/28: added ntv.R, pz.test.R, lm.jumptest.R, ntv.Rd, pz.test.Rd, lm.jumptest.Rd
- modified mpv.Rd, bns.test.Rd
\ No newline at end of file
+ modified mpv.Rd, bns.test.Rd
+2021/12/15: added ae.R, ae.cpp, ae.Rd, aeCharacteristic.Rd, aeDensity.Rd, aeExpectation.Rd, aeKurtosis.Rd, aeMarginal.Rd
+ aeMean.Rd, aeMoment.Rd, aeSd.Rd, aeSkewness.Rd, yuima.ae-class.Rd, updated NAMESPACE, updated yuima_init.c,
+ updated DESCRIPTION to add package calculus in Imports
Added: pkg/yuima/R/ae.R
===================================================================
--- pkg/yuima/R/ae.R (rev 0)
+++ pkg/yuima/R/ae.R 2021-12-15 13:53:54 UTC (rev 766)
@@ -0,0 +1,1295 @@
+#' Class for the asymptotic expansion of diffusion processes
+#'
+#' The \code{yuima.ae} class is used to describe the output of the functions \code{\link{ae}} and \code{\link{aeMarginal}}.
+#'
+#' @slot order integer. The order of the expansion.
+#' @slot var character. The state variables.
+#' @slot u.var character. The variables of the characteristic function.
+#' @slot eps.var character. The perturbation variable.
+#' @slot characteristic expression. The characteristic function.
+#' @slot density expression. The probability density function.
+#' @slot Z0 numeric. The solution to the deterministic process obtained by setting the perturbation to zero.
+#' @slot Mu numeric. The drift vector for the representation of Z1.
+#' @slot Sigma matrix. The diffusion matrix for the representation of Z1.
+#' @slot c.gamma list. The coefficients of the Hermite polynomials.
+#' @slot h.gamma list. Hermite polynomials.
+#'
+setClass("yuima.ae", slots = c(
+ order = "integer",
+ var = "character",
+ u.var = "character",
+ eps.var = "character",
+ characteristic = "expression",
+ density = "expression",
+ Z0 = "numeric",
+ Mu = "numeric",
+ Sigma = "matrix",
+ c.gamma = "list",
+ h.gamma = "list"
+))
+
+#' Constructor for yuima.ae
+#' @rdname yuima.ae-class
+setMethod("initialize", "yuima.ae", function(.Object, order, var, u.var, eps.var, characteristic, Z0, Mu, Sigma, c.gamma, verbose){
+
+ ########################################################
+ # Set density #
+ ########################################################
+ if(verbose) {
+ cat(paste0('Computing distribution...'))
+ time <- Sys.time()
+ }
+
+ tmp <- calculus::hermite(var = var, sigma = solve(Sigma), order = 3*order)
+ h.gamma <- lapply(tmp, function(x) x$f)
+ .Object at density <- sapply(0:order, function(m) {
+
+ if(m>0){
+ pdf <- sapply(1:m, function(m){
+ g <- c.gamma[[m]][sapply(c.gamma[[m]], function(x) x!=0)]
+ if(length(g)==0) return(0)
+ h <- h.gamma[names(g)]
+ p <- paste(calculus::wrap(g), calculus::wrap(h), sep = "*", collapse = " + ")
+ paste(paste0(eps.var, "^", m), calculus::wrap(p), sep = " * ")
+ })
+ pdf <- paste(pdf, collapse = " + ")
+ pdf <- paste0(1, " + ", pdf)
+ }
+ else {
+ pdf <- 1
+ }
+
+ kernel <- sprintf('%s * exp(%s)', ((2*pi)^(-length(var)/2)*(det(Sigma)^(-0.5))), (-0.5 * solve(Sigma)) %inner% (var %outer% var))
+
+ pdf <- paste0(kernel, " * (", pdf, ")")
+
+ for(i in 1:length(var)) {
+ z <- sprintf("(((%s - %s)/%s) - %s)", var[i], Z0[var[i]], eps.var, Mu[i])
+ pdf <- gsub(x = pdf, pattern = paste0("\\b",var[i],"\\b"), replacement = z)
+ pdf <- sprintf("(%s) / abs(%s)", pdf, eps.var)
+ }
+
+ return(parse(text = pdf))
+
+ })
+
+ if(verbose) {
+ cat(sprintf(' (%s sec)\n', difftime(Sys.time(), time, units = "secs")[[1]]))
+ time <- Sys.time()
+ }
+
+ ########################################################
+ # Set object #
+ ########################################################
+ .Object at order <- order
+ .Object at var <- var
+ .Object at u.var <- u.var
+ .Object at eps.var <- eps.var
+ .Object at characteristic <- characteristic
+ .Object at Z0 <- Z0
+ .Object at Mu <- Mu
+ .Object at Sigma <- Sigma
+ .Object at c.gamma <- c.gamma
+ .Object at h.gamma <- h.gamma
+
+ # return
+ return(.Object)
+
+})
+
+#' Plot method for an object of class yuima.ae
+#' @rdname yuima.ae-class
+setMethod("plot", signature(x = "yuima.ae"), function(x, grids = list(), eps = 1, order = NULL, ...){
+
+ n <- length(x at var)
+ par(mfrow = c(n,1))
+
+ for(z in x at var){
+ margin <- aeMarginal(ae = x, var = z)
+ grid <- grids[z]
+ if(is.null(grid[[1]])){
+ mu <- aeMean(ae = margin, eps = eps, order = order)[[1]]
+ sd <- aeSd(ae = margin, eps = eps, order = order)[[1]]
+ grid <- list(seq(mu-5*sd, mu+5*sd, length.out = 1000))
+ names(grid) <- z
+ }
+ dens <- do.call('aeDensity', c(grid, list(ae = margin, eps = eps, order = order)))
+ plot(x = grid[[1]], y = dens, type = 'l', xlab = z, ylab = 'Density', ...)
+ }
+
+ par(mfrow=c(1,1))
+
+})
+
+#' Asymptotic Expansion
+#'
+#' Asymptotic expansion of uni-dimensional and multi-dimensional diffusion processes.
+#'
+#' @param model an object of \code{\link{yuima-class}} or \code{\link{yuima.model-class}}.
+#' @param xinit initial value vector of state variables.
+#' @param order integer. The asymptotic expansion order. Higher orders lead to better approximations but longer computational times.
+#' @param true.parameter named list of parameters.
+#' @param sampling a \code{\link{yuima.sampling-class}} object.
+#' @param eps.var character. The perturbation variable.
+#' @param solver the solver for ordinary differential equations. One of \code{"rk4"} (more accurate) or \code{"euler"} (faster).
+#' @param verbose logical. Print on progress? Default \code{FALSE}.
+#'
+#' @return An object of \code{\link{yuima.ae-class}}
+#'
+#' @details
+#' If \code{sampling} is not provided, then \code{model} must be an object of \code{\link{yuima-class}} with non-empty \code{sampling}.
+#'
+#' if \code{eps.var} does not appear in the model specification, then it is internally added in front of the diffusion matrix to apply the asymptotic expansion scheme.
+#'
+#' @author
+#' Emanuele Guidotti <emanuele.guidotti at unine.ch>
+#'
+#' @examples
+#'
+#' # model
+#' gbm <- setModel(drift = 'mu*x', diffusion = 'sigma*x', solve.variable = 'x')
+#'
+#' # settings
+#' xinit <- 100
+#' par <- list(mu = 0.01, sigma = 0.2)
+#' sampling <- setSampling(Initial = 0, Terminal = 1, n = 1000)
+#'
+#' # asymptotic expansion
+#' approx <- ae(model = gbm, sampling = sampling, order = 4, true.parameter = par, xinit = xinit)
+#'
+#' # exact density
+#' x <- seq(50, 200, by = 0.1)
+#' exact <- dlnorm(x = x, meanlog = log(xinit)+(par$mu-0.5*par$sigma^2)*1, sdlog = par$sigma*sqrt(1))
+#'
+#' # compare
+#' plot(x, exact, type = 'l', ylab = "Density")
+#' lines(x, aeDensity(x = x, ae = approx, order = 1), col = 2)
+#' lines(x, aeDensity(x = x, ae = approx, order = 2), col = 3)
+#' lines(x, aeDensity(x = x, ae = approx, order = 3), col = 4)
+#' lines(x, aeDensity(x = x, ae = approx, order = 4), col = 5)
+#'
+#' @importFrom calculus %dot%
+#' @importFrom calculus %inner%
+#' @importFrom calculus %mx%
+#' @importFrom calculus %outer%
+#' @importFrom calculus %prod%
+#' @importFrom calculus %sum%
+#' @importFrom calculus %gradient%
+#'
+#' @export
+#'
+ae <- function(model, xinit, order = 1L, true.parameter = list(), sampling = NULL, eps.var = 'eps', solver = "rk4", verbose = FALSE){
+
+ obj.class <- class(model)
+
+ if(!is.null(sampling)){
+
+ if(obj.class=='yuima'){
+ stop('model must be of class yuima.model when sampling is provided')
+ }
+ else if(obj.class=='yuima.model'){
+ model <- setYuima(model = model, sampling = sampling)
+ }
+ else stop('model must be of class yuima or yuima.model')
+
+ }
+ else if(obj.class!='yuima'){
+ stop('model must be of class yuima when sampling is not provided')
+ }
+
+ if(!model at sampling@regular)
+ stop('ae needs regular sampling')
+
+ if(length(sampling at grid)!=1)
+ stop('ae needs a unidimensional sampling grid')
+
+ if(length(model at model@jump.coeff)>0)
+ stop('ae does not support jump processes')
+
+ if(model at model@hurst!=0.5)
+ stop('ae does not support fractional processes')
+
+ # temporary disable calculus.auto.wrap
+ calculus.auto.wrap <- options(calculus.auto.wrap = FALSE)
+ on.exit(options(calculus.auto.wrap), add = TRUE)
+
+
+ ########################################################
+ # dz #
+ ########################################################
+ dz <- function(i){
+
+ paste0('d', i, '.', AE$z)
+
+ }
+
+ dz.nu <- function(i, nu){
+
+ dz <- dz(i)
+
+ z <- sapply(1:AE$d, function(i) rep(dz[i], nu[i]), simplify = F)
+ z <- paste(unlist(z), collapse = ' * ')
+
+ if(z=='') z <- '1'
+ return(z)
+
+ }
+
+ ########################################################
+ # Z.I #
+ ########################################################
+ Z.I <- function(I, bar = FALSE) {
+
+ tmp <- NULL
+ for(i in I){
+
+ z.i <- dz(i = i)
+ if(bar) {
+ z.i <- c( z.i, ifelse(i==1, 1, 0) )
+ }
+ z.i <- array(z.i)
+
+ if(is.null(tmp)) tmp <- z.i
+ else tmp <- tmp %outer% z.i
+
+ }
+
+ return(tmp)
+
+ }
+
+ # Hermite Valued Expectation
+ HVE <- function(nu, K){
+
+ # convert to label
+ nu <- label(nu)
+
+ # for each nu'
+ H <- sapply(names(AE$c.nu[[nu]]), function(nu.prime) {
+ c(AE$c.nu[[nu]][[nu.prime]]) * as.numeric(AE$Ez.T[AE$Ez.K[[label(K)]][[nu.prime]]])
+ })
+
+ if(is.null(dim(H)))
+ H <- sum(H)
+ else
+ H <- rowSums(H)
+
+ return(H)
+ }
+
+ # Tensor Valued Expectation
+ TVE <- function(K){
+
+ K <- K + 1
+
+ nu <- AE$nu[[sum(K)]]
+
+ E <- apply(nu, 2, function(nu){
+ c <- (1i)^sum(nu) / prod(factorial(nu)) * HVE(nu = nu, K = K)
+ c <- c/prod(factorial(K))
+ paste0(AE$u,"^",nu, collapse = "*") %prod% array(calculus::wrap(c))
+ })
+
+ if(is.null(dim(E))) E <- paste(E, collapse = ' + ')
+ else E <- apply(E, 1, function(x) paste(x, collapse = ' + '))
+
+ E <- array(E, dim = rep(AE$d, length(K)))
+
+ return(E)
+ }
+
+ # Characteristic function
+ psi <- function(m){
+
+ martingale <- sprintf('exp(%s)', (calculus::wrap(1i*AE$Mu) %inner% AE$u) %sum% (-0.5 * AE$Sigma) %inner% (AE$u %outer% AE$u))
+
+ if(m>0){
+ psi <- paste(paste0(AE$eps.var, "^", (1:m)), calculus::wrap(AE$P.m[1:m]), sep = " * ", collapse = " + ")
+ psi <- paste0(1, " + ", psi)
+ }
+ else {
+ psi <- 1
+ }
+
+ psi <- paste0(martingale," * (", psi, ")")
+
+ for(i in 1:AE$d)
+ psi <- gsub(x = psi, pattern = paste0("\\b",AE$u[i],"\\b"), replacement = calculus::wrap(paste(AE$eps.var,"*",AE$u[i])))
+
+ psi <- paste0("exp((",1i,") * (",AE$u %inner% AE$Ez.T[AE$z],")) * (", psi, ")")
+
+ return(parse(text = psi))
+
+ }
+
+ # label for partitions
+ label <- function(I){
+
+ paste(I, collapse = ',')
+
+ }
+
+ ########################################################
+ # AE environment #
+ ########################################################
+ AE <- new.env()
+
+ # expansion order
+ AE$m <- as.integer(order)
+
+ # expansion variable
+ AE$eps.var <- eps.var
+
+ # model parameters
+ AE$par <- true.parameter
+ AE$par[[AE$eps.var]] <- 0
+
+ # model dimension
+ AE$d <- model at model@equation.number
+
+ # noise dimension
+ AE$r <- model at model@noise.number + 1
+
+ # solve variables
+ AE$z <- model at model@solve.variable
+
+ # extended solve variables
+ AE$z.bar <- c(AE$z, AE$eps.var)
+
+ # characteristic function variables
+ AE$u <- paste0('u', 1:AE$d)
+
+ # simulation initial value
+ AE$xinit <- xinit
+
+ # timing if verbose
+ if(verbose) time <- Sys.time()
+
+ ########################################################
+ # V: V[,1] -> V_{0}; V[,2] -> V_{1} #
+ ########################################################
+ AE$V <- NULL
+
+ # drift
+ drift <- calculus::e2c(model at model@drift)
+
+ # diffusion
+ diffusion <- unlist(lapply(model at model@diffusion, calculus::e2c))
+ diffusion <- as.array(matrix(diffusion, nrow = AE$d, ncol = AE$r-1, byrow = TRUE))
+
+ # expansion coefficient
+ is.eps.diffusion <- any(grepl(x = diffusion, pattern = paste0('\\b',AE$eps.var,'\\b')))
+ is.eps.drift <- any(grepl(x = drift, pattern = paste0('\\b',AE$eps.var,'\\b')))
+
+ if(!is.eps.diffusion){
+
+ diffusion <- AE$eps.var %prod% calculus::wrap(diffusion)
+
+ } else {
+
+ test <- parse(text = diffusion)
+
+ env <- AE$par
+ env[[AE$eps.var]] <- 0
+ for(z in AE$z)
+ env[[z]] <- runif(n = 100, min = -999, max = 999)
+
+ is.ok <- sapply(test, function(expr){
+ all(eval(expr, env)==0)
+ })
+
+ if(!all(is.ok)){
+ stop('diffusion must vanish when evaluated at epsilon = 0')
+ }
+
+ }
+
+ # building V...
+ AE$V <- array(c(drift, diffusion), c(AE$d, AE$r))
+
+ ########################################################
+ # dV #
+ ########################################################
+ AE$dV <- list()
+ if(verbose) cat('Computing dV...')
+
+ # parse v only once
+ expr <- array(parse(text = AE$V), dim = dim(AE$V))
+
+ # for j up to twice the expnasion order...
+ for(j in 1:(2*AE$m)) {
+
+ # differentiate expression
+ expr <- calculus::derivative(expr, var = AE$z.bar, deparse = FALSE)
+
+ # convert to char
+ tmp <- calculus::e2c(expr)
+
+ # break if dV vanishes
+ if(all(tmp=="0")) break
+
+ # evaluate at eps = 0
+ if(!is.eps.diffusion & !is.eps.drift){
+
+ # auto: eps * diffusion -> boost performance
+ zero <- grepl(x = tmp, pattern = paste0('\\b',AE$eps.var,'\\b'))
+ tmp[zero] <- "0"
+
+ }
+ else {
+
+ # user defined: gsub
+ tmp[] <- gsub(x = tmp, pattern = paste0('\\b',AE$eps.var,'\\b'), replacement = "0")
+
+ }
+
+ # drop white spaces (!important)
+ tmp[] <- gsub(x = tmp, pattern = ' ', replacement = '', fixed = T)
+
+ # store dV
+ AE$dV[[j]] <- tmp
+
+ }
+
+ if(verbose) {
+ cat(sprintf(' %s derivatives (%s sec)\n', length(unlist(AE$dV)), difftime(Sys.time(), time, units = "secs")[[1]]))
+ time <- Sys.time()
+ }
+
+ ########################################################
+ # dZ #
+ ########################################################
+ AE$dZ <- list()
+ if(verbose) cat('Computing dZ...')
+
+ # for j up to twice the expnasion order...
+ for(k in 1:(2*AE$m)) {
+
+ # get partitions of k
+ I.set <- calculus::partitions(n = k)
+
+ # building dZ...
+ AE$dZ[[k]] <- "0"
+
+ # for each I in I.set
+ lapply(I.set, function(I){
+
+ # length I
+ j <- length(I)
+
+ # if dV[[j]] does not vanish
+ if(j <= length(AE$dV)){
+
+ # compute U
+ U <- (factorial(sum(I))/prod(factorial(c( I, table(I) )))) %prod% calculus::wrap(AE$dV[[j]])
+
+ # isolate coefficients
+ U[] <- paste0('{',U,'}')
+
+ # add to dZ
+ AE$dZ[[k]] <- AE$dZ[[k]] %sum% ( U %dot% Z.I(I = I, bar = TRUE) )
+
+ }
+
+ # void
+ return(NULL)
+
+ })
+
+ }
+
+ if(verbose) {
+ cat(sprintf(' (%s sec)\n', difftime(Sys.time(), time, units = "secs")[[1]]))
+ time <- Sys.time()
+ }
+
+ ########################################################
+ # K.set #
+ ########################################################
+ AE$K.set <- NULL
+
+ # for n up to 4 times the expansion order...
+ for(n in 1:(4*AE$m)){
+
+ # store K.set
+ AE$K.set <- c(AE$K.set, calculus::partitions(n = n, max = 2*AE$m))
+
+ }
+
+ ########################################################
+ # Z.K #
+ ########################################################
+ AE$Z.K <- lapply(AE$K.set, function(K) Z.I(I = K))
+ names(AE$Z.K) <- lapply(AE$K.set, label)
+
+ ########################################################
+
+ if(verbose) cat('Computing Ito ODE system...')
+
+ ito <- cpp_ito(K_set = AE$K.set, dZ = AE$dZ, Z_K = AE$Z.K, d = AE$d, r = AE$r)
+ AE$ito.lhs <- ito$lhs
+ AE$ito.rhs <- ito$rhs
+ AE$ito.rhs.var <- ito$rhs.var
+
+ if(verbose) {
+ cat(sprintf(' %s equations (%s sec)\n', length(AE$ito.rhs)+AE$d, difftime(Sys.time(), time, units = "secs")[[1]]))
+ time <- Sys.time()
+ }
+
+ ########################################################
+
+ if(verbose) cat('Reducing Ito ODE system...')
+
+ AE$nu <- list()
+ AE$Ez.K <- list()
+
+ # nu indices
+ for(k in 1:(2*AE$m))
+ AE$nu[[k]] <- calculus::partitions(n = k, length = AE$d, fill = TRUE, perm = TRUE, equal = FALSE)
+
+ # find needed terms
+ for(i in 1:AE$m) {
+ for(l in 1:i){
+
+ K.set <- calculus::partitions(n = i, length = l, perm = TRUE)
+
+ lapply(K.set, function(K) {
+
+ K <- K+1
+
+ AE$Ez.K[[label(K)]] <- list()
+
+ apply(AE$nu[[sum(K)]], 2, function(nu) {
+ # store
+ AE$Ez.K[[label(K)]][[label(nu)]] <- cpp_E( dz.nu(i = 1, nu = nu) %prod% Z.I(I = K) )
+ # void
+ return(NULL)
+ })
+
+ # void
+ return(NULL)
+ })
+ }
+ }
+
+ ########################################################
+
+ AE$Ez.T <- list()
+ AE$Ez <- unique(unlist(AE$Ez.K))
+
+ # drop duplicated equations
+ idx <- which(!duplicated(AE$ito.lhs))
+ AE$ito.lhs <- AE$ito.lhs[idx]
+ AE$ito.rhs <- AE$ito.rhs[idx]
+ AE$ito.rhs.var <- AE$ito.rhs.var[idx]
+
+ while(TRUE){
+ # needed equation id
+ idx <- which(AE$ito.lhs %in% AE$Ez)
+ # additional needed var
+ add <- unique(unlist(AE$ito.rhs.var[idx]))
+ add <- add[!(add %in% AE$Ez)]
+ # break or store
+ if(length(add)==0) break
+ AE$Ez <- c(AE$Ez, add)
+ }
+
+ # drop equations not needed
+ idx <- which(AE$ito.lhs %in% AE$Ez)
+ AE$ito.lhs <- AE$ito.lhs[idx]
+ AE$ito.rhs <- AE$ito.rhs[idx]
+ AE$ito.rhs.var <- AE$ito.rhs.var[idx]
+
+ # plug zero if empty equation
+ while(TRUE){
+
+ idx <- which(AE$ito.rhs=='')
+ if(length(idx)==0) break
+
+ zero <- AE$ito.lhs[idx]
+ AE$Ez.T[zero] <- 0
+ AE$ito.lhs <- AE$ito.lhs[-idx]
+ AE$ito.rhs <- AE$ito.rhs[-idx]
+
+ pattern <- gsub(zero, pattern = '.', replacement = '\\.', fixed = T)
+ pattern <- gsub(pattern, pattern = '_', replacement = '\\_', fixed = T)
+ pattern <- paste0(pattern, collapse = '|')
+ pattern <- paste0(' \\* (',pattern,')\\b')
+
+ AE$ito.rhs <- unlist(lapply(strsplit(AE$ito.rhs, split = ' + ', fixed = T), function(x){
+ is.zero <- grepl(x = x, pattern = pattern)
+ x <- x[!is.zero]
+ if(length(x)==0) return("")
+ else return(paste(x, collapse = ' + '))
+ }))
+
+ }
+
+ if(verbose) {
+ cat(sprintf(' %s equations (%s sec)\n', length(AE$ito.rhs)+AE$d, difftime(Sys.time(), time, units = "secs")[[1]]))
+ time <- Sys.time()
+ }
+
+ ########################################################
+
+ if(verbose) cat('Solving Ito ODE system...')
+
+ # Ito ODE
+ lhs <- c(AE$z, AE$ito.lhs)
+ rhs <- c(AE$V[,1], AE$ito.rhs)
+
+ # Initial values
+ xinit <- c(rep(AE$xinit, length.out = AE$d), rep(0, length(lhs)-AE$d))
+ names(xinit) <- lhs
+
+ # Solve
+ AE$Ez.T <- calculus::ode(f = rhs, var = xinit, times = sampling at grid[[1]], params = AE$par, timevar = model at model@time.variable, drop = TRUE, method = solver)
+
+ if(verbose) {
+ cat(sprintf(' (%s sec)\n', difftime(Sys.time(), time, units = "secs")[[1]]))
+ time <- Sys.time()
+ }
+
+ ########################################################
+
+ if(verbose) cat('Computing Sigma matrix...')
+
+ y.lhs <- array(paste('y', gsub(AE$z %outer% AE$z, pattern = " * ", replacement = "_", fixed = T), sep = '.'), dim = rep(AE$d, 2))
+ y.rhs <- calculus::wrap(AE$V[,1] %gradient% AE$z) %mx% y.lhs
+ y.0 <- diag(AE$d)
+
+ dim(y.lhs) <- NULL
+ dim(y.rhs) <- NULL
+ dim(y.0) <- NULL
+
+ y.lhs <- c(AE$z, y.lhs)
+ y.rhs <- c(AE$V[,1], y.rhs)
+
+ y.0 <- c(rep(AE$xinit, length.out = AE$d), y.0)
+ names(y.0) <- y.lhs
+
+ times <- sampling at grid[[1]]
+ n.times <- length(times)
+
+ x <- calculus::ode(f = y.rhs, var = y.0, times = times, params = AE$par, timevar = model at model@time.variable, method = solver)
+ y <- x[, -c(1:AE$d), drop = FALSE]
+ z <- as.data.frame(x[, c(1:AE$d), drop = FALSE])
+ if(!is.null(model at model@time.variable))
+ z[[model at model@time.variable]] <- times
+
+ y.s <- lapply(1:n.times, function(i) matrix(y[i,], nrow = AE$d))
+ y_inv.s <- lapply(y.s, solve)
+
+ dV <- AE$V %gradient% AE$eps.var
+ dV.s <- calculus::evaluate(dV, as.data.frame(c(z, AE$par)))
+ dV.s <- lapply(1:n.times, function(i) array(dV.s[i,], dim = dim(dV)))
+ if(length(dV.s)==1)
+ dV.s = rep(dV.s, n.times)
+
+ # Mu
+ if(all(dV[,1,]=="0")){
+
+ AE$Mu <- array(0, dim = AE$d)
+
+ }
+ else {
+
+ Mu <- sapply(1:n.times, function(i){
+ y_inv.s[[i]] %*% dV.s[[i]][,1,]
+ })
+
+ if(is.null(dim(Mu))) {
+ n <- length(Mu)
+ Mu <- ( sum(Mu[c(-1,-n)]) + sum(Mu[c(1,n)])/2 ) * sampling at delta
+ }
+ else {
+ n <- ncol(Mu)
+ Mu <- ( rowSums(Mu[,c(-1,-n)]) + rowSums(Mu[,c(1,n)])/2 ) * sampling at delta
+ }
+
+ AE$Mu <- array(y.s[[n.times]] %*% Mu)
+
+ }
+
+ # Sigma
+ S <- sapply(1:n.times, function(i){
+ a <- y.s[[n.times]] %*% y_inv.s[[i]] %*% dV.s[[i]][,-1,]
+ a %*% t(a)
+ })
+
+ if(is.null(dim(S))) {
+ n <- length(S)
+ S <- ( sum(S[c(-1,-n)]) + sum(S[c(1,n)])/2 ) * sampling at delta
+ }
+ else {
+ n <- ncol(S)
+ S <- ( rowSums(S[,c(-1,-n)]) + rowSums(S[,c(1,n)])/2 ) * sampling at delta
+ }
+
+ AE$Sigma <- array(S, dim = c(AE$d,AE$d))
+
+ if(verbose) {
+ cat(sprintf(' (%s sec)\n', difftime(Sys.time(), time, units = "secs")[[1]]))
+ time <- Sys.time()
+ }
+
+ ########################################################
+
+ # Hermite coefficients
+ if(verbose) cat(paste0('Computing Hermite...'))
+
+ tmp <- calculus::hermite(var = AE$z, sigma = AE$Sigma, order = 2*AE$m,
+ transform = solve(AE$Sigma) %dot% calculus::wrap(AE$z %sum% -AE$Mu))
+
+ AE$c.nu <- lapply(tmp, function(x) {
+ coef <- as.list(x$terms$coef)
+ names(coef) <- rownames(x$terms)
+ return(coef)
+ })
+ AE$h.nu <- lapply(tmp, function(x) x$f)
+
+ if(verbose) {
+ cat(sprintf(' (%s sec)\n', difftime(Sys.time(), time, units = "secs")[[1]]))
+ time <- Sys.time()
+ }
+
+ ########################################################
+
+ AE$ul <- list(array(AE$u))
+ if(AE$m > 1) for(l in 2:AE$m){
+ AE$ul[[l]] <- AE$ul[[l-1]] %outer% AE$u
+ }
+ AE$ul <- lapply(AE$ul, function(ul){
+ array(unlist(lapply(strsplit(ul, split = " * ", fixed = T), function(u) {
+ u <- table(u)
+ paste0(names(u),"^",u, collapse = "*")
+ })), dim = dim(ul))
+ })
+
+ ########################################################
+
+ if(verbose) cat('Computing characteristic function...')
+
+ AE$psi <- list()
+ for(m in 1:AE$m) {
+
+ AE$psi[[m]] <- list()
+
+ for(l in 1:m){
+
+ K.set <- calculus::partitions(n = m, length = l, perm = TRUE)
+
+ psi.m.l <- unlist(lapply(K.set, function(K){
+ sprintf("(%s) * (%s)", (1i)^l, calculus::wrap(TVE(K = K)) %inner% AE$ul[[l]])
+ }))
+
+ expr <- sprintf('1/%s * (%s)', factorial(l), paste(psi.m.l, collapse = ' + '))
+ AE$psi[[m]][[l]] <- calculus::taylor(expr, var = AE$u, order = m+2*l)$f
+
+ }
+
+ }
+
+ AE$P.m = sapply(AE$psi, function(p.m.l) paste(p.m.l, collapse = " + "))
+ AE$c.gamma <- lapply(1:AE$m, function(m) {
+ p <- calculus::taylor(AE$P.m[m], var = AE$u, order = 3*m)
+ coef <- Re(p$terms$coef/(1i)^p$terms$degree)
+ coef <- as.list(coef)
+ names(coef) <- rownames(p$terms)
+ return(coef)
+ })
+
+ AE$PSI <- sapply(0:AE$m, psi)
+
+ if(verbose) {
+ cat(sprintf(' (%s sec)\n', difftime(Sys.time(), time, units = "secs")[[1]]))
+ time <- Sys.time()
+ }
+
+ ########################################################
+
+ return(new(
+ "yuima.ae",
+ order = AE$m,
+ var = AE$z,
+ u.var = AE$u,
+ eps.var = AE$eps.var,
+ characteristic = AE$PSI,
+ Z0 = unlist(AE$Ez.T[AE$z]),
+ Mu = as.numeric(AE$Mu),
+ Sigma = AE$Sigma,
+ c.gamma = AE$c.gamma,
+ verbose = verbose
+ ))
+
+}
+
+#' Asymptotic Expansion - Density
+#'
+#' @param ... named argument, data.frame, list, or environment specifying the grid to evaluate the density. See examples.
+#' @param ae an object of class \code{\link{yuima.ae-class}}.
+#' @param eps numeric. The intensity of the perturbation.
+#' @param order integer. The expansion order. If \code{NULL} (default), it uses the maximum order used in \code{ae}.
+#'
+#' @return Probability density function evaluated on the given grid.
+#'
+#' @examples
+#'
+#' # model
+#' gbm <- setModel(drift = 'mu*x', diffusion = 'sigma*x', solve.variable = 'x')
+#'
+#' # settings
+#' xinit <- 100
+#' par <- list(mu = 0.01, sigma = 0.2)
+#' sampling <- setSampling(Initial = 0, Terminal = 1, n = 1000)
+#'
+#' # asymptotic expansion
+#' approx <- ae(model = gbm, sampling = sampling, order = 4, true.parameter = par, xinit = xinit)
+#'
+#' # The following are all equivalent methods to specify the grid via ....
+#' # Notice that the character 'x' corresponds to the solve.variable of the yuima model.
+#'
+#' # 1) named argument
+#' x <- seq(50, 200, by = 0.1)
+#' density <- aeDensity(x = x, ae = approx, order = 4)
+#' # 2) data frame
+#' df <- data.frame(x = seq(50, 200, by = 0.1))
+#' density <- aeDensity(df, ae = approx, order = 4)
+#' # 3) environment
+#' env <- new.env()
+#' env$x <- seq(50, 200, by = 0.1)
+#' density <- aeDensity(env, ae = approx, order = 4)
+#' # 4) list
+#' lst <- list(x = seq(50, 200, by = 0.1))
+#' density <- aeDensity(lst, ae = approx, order = 4)
+#'
+#' # exact density
+#' exact <- dlnorm(x = x, meanlog = log(xinit)+(par$mu-0.5*par$sigma^2)*1, sdlog = par$sigma*sqrt(1))
+#'
+#' # compare
+#' plot(x = exact, y = density, xlab = "Exact", ylab = "Approximated")
+#'
+#' @export
+#'
+aeDensity <- function(..., ae, eps = 1, order = NULL){
+
+ if(is.null(order))
+ order <- ae at order
+
+ pdf <- ae at density[order+1]
+
+ env <- list(...)
+ if(length(env)==1)
+ if(is.list(env[[1]]) | is.environment(env[[1]]))
+ env <- env[[1]]
+
+ env[[ae at eps.var]] <- eps
+
+ return(eval(expr = pdf, envir = env))
+
+}
+
+#' Asymptotic Expansion - Marginals
+#'
+#' @param ae an object of class \code{\link{yuima.ae-class}}.
+#' @param var variables of the marginal distribution to compute.
+#'
+#' @return An object of \code{\link{yuima.ae-class}}
+#'
+#' @examples
+#'
+#' # multidimensional model
+#' gbm <- setModel(drift = c('mu*x1','mu*x2'),
+#' diffusion = matrix(c('sigma1*x1',0,0,'sigma2*x2'), nrow = 2),
+#' solve.variable = c('x1','x2'))
+#'
+#' # settings
+#' xinit <- c(100, 100)
+#' par <- list(mu = 0.01, sigma1 = 0.2, sigma2 = 0.1)
+#' sampling <- setSampling(Initial = 0, Terminal = 1, n = 1000)
+#'
+#' # asymptotic expansion
+#' approx <- ae(model = gbm, sampling = sampling, order = 3, true.parameter = par, xinit = xinit)
+#'
+#' # extract marginals
+#' margin1 <- aeMarginal(ae = approx, var = "x1")
+#' margin2 <- aeMarginal(ae = approx, var = "x2")
+#'
+#' # compare with exact solution for marginal 1
+#' x1 <- seq(50, 200, by = 0.1)
+#' exact <- dlnorm(x = x1, meanlog = log(xinit[1])+(par$mu-0.5*par$sigma1^2), sdlog = par$sigma1)
+#' plot(x1, exact, type = 'p', ylab = "Density")
+#' lines(x1, aeDensity(x1 = x1, ae = margin1, order = 3), col = 2)
+#'
+#' # compare with exact solution for marginal 2
+#' x2 <- seq(50, 200, by = 0.1)
+#' exact <- dlnorm(x = x2, meanlog = log(xinit[2])+(par$mu-0.5*par$sigma2^2), sdlog = par$sigma2)
+#' plot(x2, exact, type = 'p', ylab = "Density")
+#' lines(x2, aeDensity(x2 = x2, ae = margin2, order = 3), col = 2)
+#'
+#' @export
+#'
+aeMarginal <- function(ae, var){
+
+ # init
+ keep <- ae at var %in% var
+ if(sum(!keep)==0) return(ae)
+
+ # vanish marginal coefficients
+ c.gamma <- ae at c.gamma
+ for(i in 1:length(c.gamma)) for(j in 1:length(c.gamma[[i]])){
+ o <- as.numeric(unlist(strsplit(x = names(c.gamma[[i]][j]), split = ',')))
+ if(any(o[!keep]>0)){
+ c.gamma[[i]][[j]] <- 0
+ }
+ else {
+ names(c.gamma[[i]])[j] <- paste(o[keep], collapse = ',')
+ }
+ }
+
+ # characteristic
+ characteristic <- sapply(calculus::e2c(ae at characteristic), function(psi) {
+ for(u in ae at u.var[!keep])
+ psi <- gsub(x = psi, pattern = paste0("\\b",u,"\\b"), replacement = 0)
+ return(calculus::c2e(psi))
+ })
+
+ # return
+ return(new(
+ "yuima.ae",
+ order = ae at order,
+ var = ae at var[keep],
+ u.var = ae at u.var[keep],
+ eps.var = ae at eps.var,
+ characteristic = characteristic,
+ Z0 = ae at Z0[keep],
+ Mu = ae at Mu[keep],
+ Sigma = ae at Sigma[keep, keep, drop = FALSE],
+ c.gamma = c.gamma,
+ verbose = FALSE
+ ))
+
+}
+
+#' Asymptotic Expansion - Functionals
+#'
+#' Compute the expected value of functionals.
+#'
+#' @param f character. The functional.
+#' @param bounds named list of integration bounds in the form \code{list(x = c(xmin, xmax), y = c(ymin, ymax), ...)}
+#' @param ae an object of class \code{\link{yuima.ae-class}}.
+#' @param eps numeric. The intensity of the perturbation.
+#' @param order integer. The expansion order. If \code{NULL} (default), it uses the maximum order used in \code{ae}.
+#' @param ... additional arguments passed to \code{\link[cubature]{cubintegrate}}.
+#'
+#' @return return value of \code{\link[cubature]{cubintegrate}}. The expectation of the functional provided.
+#'
+#' @examples
+#'
+#' # model
+#' gbm <- setModel(drift = 'mu*x', diffusion = 'sigma*x', solve.variable = 'x')
+#'
+#' # settings
+#' xinit <- 100
+#' par <- list(mu = 0.01, sigma = 0.2)
+#' sampling <- setSampling(Initial = 0, Terminal = 1, n = 1000)
+#'
+#' # asymptotic expansion
+#' approx <- ae(model = gbm, sampling = sampling, order = 4, true.parameter = par, xinit = xinit)
+#'
+#' # compute the mean via integration
+#' aeExpectation(f = 'x', bounds = list(x = c(0,1000)), ae = approx)
+#'
+#' # compare with the mean computed by differentiation of the characteristic function
+#' aeMean(approx)
+#'
+#' @export
+#'
+aeExpectation <- function(f, bounds, ae, eps = 1, order = NULL, ...){
+
+ f <- calculus::c2e(f)
+ var <- names(bounds)
+ ae <- aeMarginal(ae = ae, var = var)
+ lower <- sapply(bounds, function(x) x[1])
+ upper <- sapply(bounds, function(x) x[2])
+
+ args <- list(...)
+ args$f <- function(x) {
+ names(x) <- var
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 766
More information about the Yuima-commits
mailing list