estimate-mixedDiffusion-method: Estimation for hierarchical (mixed) diffusion model

In SimoneHermann/BaPreStoPro: Bayesian Prediction of Stochastic Processes

Description

Bayesian estimation of a model dY_t = b(φ_j,t,Y_t)dt + γ \widetilde{s}(t,Y_t)dW_t, φ_j\sim N(μ, Ω), Y_{t_0}=y_0(φ, t_0).

Usage

## S4 method for signature 'mixedDiffusion'
estimate(model.class, t, data, nMCMC, propSd,
  adapt = TRUE, proposal = c("normal", "lognormal"))

Arguments

model.class	class of the hierarchical diffusion model including all required information, see mixedDiffusion-class
t	list or vector of time points
data	list or matrix of observation variables
nMCMC	length of Markov chain
propSd	vector of proposal variances for φ
adapt	if TRUE (default), proposal variance is adapted
proposal	proposal density: "normal" (default) or "lognormal" (for positive parameters)

References

Hermann, S., K. Ickstadt and C. H. Mueller (2016). Bayesian Prediction of Crack Growth Based on a Hierarchical Diffusion Model. Applied Stochastic Models in Business and Industry, DOI: 10.1002/asmb.2175.

Examples

mu <- 2; Omega <- 0.4; phi <- matrix(rnorm(21, mu, sqrt(Omega)))
model <- set.to.class("mixedDiffusion",
             parameter = list(phi = phi, mu = mu, Omega = Omega, gamma2 = 0.1),
             b.fun = function(phi, t, x) phi*x, sT.fun = function(t, x) x)
t <- seq(0, 1, by = 0.01)
data <- simulate(model, t = t, plot.series = TRUE)
est <- estimate(model, t, data[1:20,], 100)  # nMCMC should be much larger
plot(est)

# OU
b.fun <- function(phi, t, y) phi[1]-phi[2]*y; y0.fun <- function(phi, t) phi[3]
mu <- c(10, 5, 0.5); Omega <- c(0.9, 0.01, 0.01)
phi <- sapply(1:3, function(i) rnorm(21, mu[i], sqrt(Omega[i])))
model <- set.to.class("mixedDiffusion",
               parameter = list(phi = phi, mu = mu, Omega = Omega, gamma2 = 0.1),
               y0.fun = y0.fun, b.fun = b.fun, sT.fun = function(t, x) 1)
t <- seq(0, 1, by = 0.01)
data <- simulate(model, t = t, plot.series = TRUE)
est <- estimate(model, t, data[1:20,], 100)  # nMCMC should be much larger
plot(est)

##
t.list <- list()
for(i in 1:20) t.list[[i]] <- t
t.list[[21]] <- t[1:50]
data.list <- list()
for(i in 1:20) data.list[[i]] <- data[i,]
data.list[[21]] <- data[21, 1:50]
est <- estimate(model, t.list, data.list, 100)
pred <- predict(est, t = t[50:101], which.series = "current", ind.pred = 21,
   b.fun.mat = function(phi, t, y) phi[,1]-phi[,2]*y)

SimoneHermann/BaPreStoPro documentation built on May 9, 2019, 1:46 p.m.