BaPreStoPro (version 0.1)

predict,est.mixedDiffusion-method: Prediction for a hierarchical (mixed) diffusion process model

Description

Bayesian prediction of a stochastic process model $dY_t = b(\phi_j,t,Y_t)dt + \gamma \widetilde{s}(t,Y_t)dW_t, \phi_j~N(\mu, \Omega)$.

Usage

"predict"(object, t, Euler.interval = FALSE, level = 0.05, burnIn, thinning, b.fun.mat, which.series = c("new", "current"), y.start, ind.pred, M2pred = 10, cand.length = 1000, pred.alg = c("Distribution", "Trajectory", "simpleTrajectory", "simpleBayesTrajectory"), sample.length, grid, plot.prediction = TRUE)

Arguments

object

class object of MCMC samples: "est.mixedDiffusion", created with method estimate,mixedDiffusion-method

t

vector of time points to make predictions for

Euler.interval

if TRUE: simple prediction intervals with Euler are made (in one step each)

level

level of the prediction intervals

burnIn

burn-in period

thinning

thinning rate

b.fun.mat

matrix-wise definition of drift function (makes it faster)

which.series

which series to be predicted, new one ("new") or further development of current one ("current")

y.start

optional, if missing, first (which.series = "new") or last observation variable ("current") is taken

ind.pred

index of series to be predicted, optional, if which.series = "current" and ind.pred missing, the last series is taken

M2pred

optional, if current series to be predicted and t missing, M2pred variables will be predicted with the observation time distances

cand.length

length of candidate samples (if method = "vector")

pred.alg

prediction algorithm, "Distribution", "Trajectory", "simpleTrajectory" or "simpleBayesTrajectory"

sample.length

number of samples to be drawn, default is the number of posterior samples

grid

fineness degree of sampling approximation

plot.prediction

if TRUE, prediction intervals are plotted

References

Hermann, S. (2016a). BaPreStoPro: an R Package for Bayesian Prediction of Stochastic Processes. SFB 823 discussion paper 28/16.

Hermann, S. (2016b). Bayesian Prediction for Stochastic Processes based on the Euler Approximation Scheme. SFB 823 discussion paper 27/16.

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)
est_mixdiff <- estimate(model, t, data[1:20,], 100) # nMCMC should be much larger
plot(est_mixdiff)
## Not run: 
# pred_mixdiff <- predict(est_mixdiff, b.fun.mat = function(phi, t, y) phi[,1]*y)
# lines(t, data[21,], lwd = 2)
# mean(apply(pred_mixdiff$Y, 2, quantile, 0.025) <= data[21, ] &
# apply(pred_mixdiff$Y, 2, quantile, 0.975) >= data[21, ])
# mean(sapply(1:20, function(i){
#    mean(apply(pred_mixdiff$Y, 2, quantile, 0.025) <= data[i, ] &
#    apply(pred_mixdiff$Y, 2, quantile, 0.975) >= data[i, ])}))
# pred_mixdiff2 <- predict(est_mixdiff, b.fun.mat = function(phi, t, y) phi[,1]*y,
#      which.series = "current")
# pred_mixdiff3 <- predict(est_mixdiff, b.fun.mat = function(phi, t, y) phi[,1]*y,
#      which.series = "current", y.start = data[20, 51], t = t[51:101])
# ## End(Not run)
pred_mixdiff <- predict(est_mixdiff, Euler.interval = TRUE,
     b.fun.mat = function(phi, t, y) phi[,1]*y); lines(t, data[21,], lwd = 2)
     # one step Euler approximation
pred_mixdiff <- predict(est_mixdiff, pred.alg = "simpleTrajectory",
                        sample.length = 100)
for(i in 1:100) lines(t, pred_mixdiff$Y[i,], col = "grey")
pred_mixdiff <- predict(est_mixdiff, pred.alg = "simpleBayesTrajectory")

# OU
## Not run: 
# b.fun <- function(phi, t, y) phi[1]-phi[2]*y; y0.fun <- function(phi, t) phi[3]
# mu <- c(10, 1, 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)
# est <- estimate(model, t, data[1:20,], 2000)
# plot(est)
# pred <- predict(est, t = seq(0, 1, length = 21),
#    b.fun.mat = function(phi, t, y) phi[,1]-phi[,2]*y)
# lines(t, data[21,], lwd = 2)
# mean(apply(pred$Y, 2, quantile, 0.025) <= data[21, seq(1, length(t), length = 21)] &
#      apply(pred$Y, 2, quantile, 0.975) >= data[21, seq(1, length(t), length = 21)])
# mean(sapply(1:20, function(i){
#     mean(apply(pred$Y, 2, quantile, 0.025) <= data[i, seq(1, length(t), length = 21)] &
#     apply(pred$Y, 2, quantile, 0.975) >= data[i, seq(1, length(t), length = 21)])}))
# ## End(Not run)