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

In 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)

Example output

10 of 51 predictions are calculated
20 of 51 predictions are calculated
30 of 51 predictions are calculated
40 of 51 predictions are calculated
50 of 51 predictions are calculated

BaPreStoPro documentation built on May 2, 2019, 3:34 p.m.