BaPreStoPro-package: Bayesian Prediction of Stochastic Processes

In SimoneHermann/BaPreStoPro: Bayesian Prediction of Stochastic Processes

Description

This package contains simulate, estimate and predict methods for non-homogeneous Poisson processes (NHPP), jump diffusions, (mixed) diffusions, hidden (mixed) diffusion models, regression model including a NHPP, and (mixed) regression models for comparison.

Details

Package:	BaPreStoPro
Type:	Package
Version:	1.0
Date:	2016-06-07
License:	GLP-2, GLP-3

Each of the models has its specific class, "jumpDiffusion", "Merton", "Diffusion", "mixedDiffusion", "hiddenDiffusion", "hiddenmixedDiffusion", "jumpRegression", "NHPP", "Regression", "mixedRegression", created with the function set.to.class. For each of the model classes, a method simulate and estimate are provided. The output of method estimate is a new class object with the prefix "est.". For the estimation classes, methods plot and predict are available. An overview of the package can be found in Hermann (2016a) and theoretical details to the prediction procedures in Hermann (2016b).

jumpDiffusion

The model class "jumpDiffusion" contains all information about model defined by the stochastic differential equation (SDE) dY_t = b(φ,t,Y_t)dt + s(γ^2,t,Y_t)dW_t + h(θ,t,Y_t)dN_t with N_t\sim Pois(Λ(t, ξ)) a non-homogeneous Poisson process and W_t a Brownian motion. The SDE is approximated with the Euler Maruyama approximation, which leads, dependent on the Poisson process variables, to a normal likelihood. For more information how to build the class, see examples in jumpDiffusion-class and estimate,jumpDiffusion-method.

Estimation is done by a Metropolis-within-Gibbs sampler. For each of the parameters, a Metropolis-Hastings (MH) step is made, where the proposal density can be chosen between normal and lognormal. A proposal standard deviation can be chosen, which, if desired, is adapted after every 50 iterations, see Rosenthal (2011).

In the case of unobserved variables of the Poisson process, one step of the Gibbs sampler is filtering of the unobserved jump process. Details can be found in Hermann and Ruggeri (2016) or Hermann (2016a).

Merton

Specific choices of functions b, s and h lead to an explicit solution of the process Y_t = y_0 \exp( φ t - γ^2/2 t+γ W_t + \log(1+θ) N_t). This model is well-known in the literature as Merton model. There are conjugate prior distributions available for φ, γ^2 and \log(1+θ): φ and \log(1+θ) are assumed to have normal prior distribution with parameters m.phi (mean) and v.phi (variance), m.thetaT (mean) and v.thetaT (variance) respectively. γ^2 is assumed to have an inverse gamma distribution with parameters alpha.gamma and beta.gamma. An example how to build the model class can be found in Merton-class and estimate,Merton-method.

Estimation is similar to the jump diffusion process based on the Euler approximation. The difference is that for the parameters with conjugate priors, no MH step is necessary and drawing from the full conditional posterior is possible.

Diffusion

The special case of h(θ, t, y)=0 leads to a general diffusion process dY_t = b(φ,t,Y_t)dt + s(γ^2,t,Y_t)dW_t. We here restrict to the special case of s(γ^2, t, y)=√{γ^2}\widetilde{s}(t, y), because a conjugate prior for γ^2 is available (the inverse gamma with parameters alpha.gamma and beta.gamma) in this case. For φ, a normal prior with parameters m.phi (mean) and v.phi (variance) is assumed. An example can be found in Diffusion-class and estimate,Diffusion-method.

A Gibbs sampler with an MH step for φ and one step drawing from the full conditional of γ^2 is implemented. For φ, the proposal density can be chosen, "normal" or "lognormal" and the proposal standard deviation, with the option to adapt, as well.

mixedDiffusion

The diffusion process is extended to a hierarchical model with φ as random effect with normal mixture distribution with mean μ and variance Ω (diagonal matrix). For μ, a normal prior (parameters: m.mu and v.mu) is conjugate. For each diagonal element of Ω, an inverse gamma prior (parameters: alpha.omega and beta.omega) is conjugate. Further information can be found in mixedDiffusion-class and estimate,mixedDiffusion-method.

A Gibbs sampler with an MH step for each random effect and each one step drawing from the full conditionals of μ, Ω and γ^2 is implemented.

hiddenDiffusion

The same model as the diffusion above is taken, but with an added error: Z_i = Y_{t_i} + ε_i, dY_t = b(φ,t,Y_t)dt + γ \widetilde{s}(t,Y_t)dW_t, ε_i\sim N(0,σ^2), Y_{t_0}=y_0(φ, t_0). The inverse gamma prior for σ^2 is conjugate and is, therefore, implemented with parameters alpha.sigma and beta.sigma. The diffusion process is a latent variable and has also to be estimated. A conditional sequential Monte Carlo (SMC) approach is implemented. For further details see Andrieu et al. (2010). Examples can be found in hiddenDiffusion-class and estimate,hiddenDiffusion-method

hiddenmixedDiffusion

Here, the hidden diffusion model is extended to a hierarchical model with the random effect φ, similar to the mixed diffusion model. In the Gibbs sampler, one step is filtering the unobserved diffusion process for each observed series. Based on this estimation, the random effects are estimated with an MH step. With conjugate priors, full conditionals for μ, Ω, γ^2 and σ^2 are available. Examples can be found in hiddenmixedDiffusion-class and estimate,hiddenmixedDiffusion-method.

jumpRegression

We here consider the regression model: y_i = f(t_i, N_{t_i}, θ) + ε_i dependent on the Poisson process N_t\sim Pois(Λ(t, ξ)) with ε_i\sim N(0,γ^2\widetilde{s}(t)).

Here, for the case of missing observations for the Poisson process variable, a filtering procedure based on the conditional SMC is implemented. But this is still work in progress.

See examples in jumpRegression-class and estimate,jumpRegression-method.

NHPP

Some models base on the non-homogeneous Poisson process. Here, only the NHPP itself is considered. A simple MH algorithm is implemented. A proposal density can be chosen, normal or lognormal, and the proposal standard deviation as well, which is adapted, if desired.

See examples in NHPP-class and estimate,NHPP-method.

Regression

For the case of a comparison of regression and diffusion model, as made, for example, in Hermann et al. (2016), estimation and prediction is also made for regression models. All notations are analogously to the diffusion model. See examples in Regression-class and estimate,Regression-method.

mixedRegression

Analogous hierarchical regression model to the mixed diffusion model above. See examples in mixedRegression-class and estimate,mixedRegression-method.

Author(s)

Simone Hermann <[e-mail]>

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.

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.

Hermann, S. and F. Ruggeri (2016). Modelling Wear in Cylinder Liners. SFB 823 discussion paper 06/16.

Hermann, S., K. Ickstadt, and C. H. Mueller (2015). Bayesian Prediction for a Jump Diffusion Process with Application to Crack Growth in Fatigue Experiments. SFB 823 discussion paper 30/15.

Heeke, G., S. Hermann, R. Maurer, K. Ickstadt, and C. H. Mueller (2015). Stochastic Modeling and Statistical Analysis of Fatigue Tests on Prestressed Concrete Beams under Cyclic Loadings. SFB 823 discussion paper 25/15.

Monte Carlo methods:

Robert, C. P. and G. Casella (2004). Monte Carlo Statistical Methods. Springer, New York.

Adaptive MCMC:

Rosenthal, J. S. (2011). Optimal Proposal Distributions and Adaptive MCMC. In: Handbook of Markov Chain Monte Carlo, pp. 93-112.

particel Gibbs / SMC:

Andrieu, C., A. Doucet and R. Holenstein (2010). Particle Markov Chain Monte Carlo Methods. Journal of the Royal Statistical Society B 72, pp. 269-342.

Examples

model <- set.to.class("Diffusion", parameter = list(phi = 0.5, gamma2 = 0.01))
t <- seq(0, 1, by = 0.1)
data <- simulate(model, t = t, y0 = 0.5, plot.series = TRUE)
est <- estimate(model, t, data, 10) # better: 10000
plot(est)
pred <- predict(est)

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