BaPreStoPro (version 0.1)

estimate,jumpDiffusion-method: Estimation for jump diffusion process

Description

Bayesian estimation of a stochastic process $dY_t = b(\phi,t,Y_t)dt + s(\gamma^2,t,Y_t)dW_t + h(\theta,t,Y_t)dN_t$.

Usage

"estimate"(model.class, t, data, nMCMC, propSd, adapt = TRUE, proposal = c("normal", "lognormal"), it.xi = 5)

Arguments

model.class

class of the jump diffusion model including all required information, see jumpDiffusion-class

t

vector of time points

data

vector of observation variables

nMCMC

length of Markov chain

propSd

vector of proposal variances for $(\phi, \theta, \gamma^2, \xi)$

adapt

if TRUE (default), proposal variance is adapted

proposal

proposal density for phi, theta: "normal" (default) or "lognormal" (for positive parameters), see description below

it.xi

number of iterations for MH step for $\xi$ inside the Gibbs sampler

Proposal densities

For $\gamma^2$, always the lognormal density is taken, since the parameter is always positive. For $\theta$ and $\phi$, there is the possibility to choose "normal" or "lognormal" (for both together). The proposal density for $\xi$ depends on the starting value of $\xi$. If all components are positive, the proposal density is lognormal, and normal otherwise.

Examples

# non-informative
model <- set.to.class("jumpDiffusion", Lambda = function(t, xi) (t/xi[2])^xi[1],
               parameter = list(theta = 0.1, phi = 0.05, gamma2 = 0.1, xi = c(3, 1/4)))
t <- seq(0, 1, by = 0.01)
data <- simulate(model, t = t, y0 = 0.5, plot.series = TRUE)
est <- estimate(model, t, data, 1000)
plot(est)

# informative
model <- set.to.class("jumpDiffusion", Lambda = function(t, xi) (t/xi[2])^xi[1],
   parameter = list(theta = 0.1, phi = 0.05, gamma2 = 0.1, xi = c(3, 1/4)),
   priorDensity = list(phi = function(phi) dnorm(phi, 0.05, 0.01),
                       theta = function(theta) dgamma(1/theta, 10, 0.1*9),
                       gamma2 = function(gamma2) dgamma(1/gamma2, 10, 0.1*9),
                       xi = function(xi) dnorm(xi, c(3, 1/4), c(1,1))))
t <- seq(0, 1, by = 0.01)
data <- simulate(model, t = t, y0 = 0.5, plot.series = TRUE)
est <- estimate(model, t, data, 1000)
plot(est)

## Not run: 
# est_hidden <- estimate(model, t, data$Y, 1000)
# plot(est_hidden)
# ## End(Not run)