msde.fit: Estimation Of The Random Effects In Mixed Stochastic...

In MsdeParEst: Parametric Estimation in Mixed-Effects Stochastic Differential Equations

Description

Parametric estimation of the joint density of the random effects in the mixed SDE

dX_j(t)= (α_j- β_j X_j(t))dt + σ_j \ a(X_j(t)) dW_j(t),

j=1,…,M, where the (W_j(t)) are independant Wiener processes and the (X_j(t)) are observed without noise. There can be random effects either in the drift (α_j,β_j) or in the diffusion coefficient σ_j or both (α_j,β_j,σ_j).

Usage

msde.fit(times, X, model = c("OU", "CIR"), drift.random = c(1, 2),
  drift.fixed = NULL, diffusion.random = 0, diffusion.fixed = NULL,
  nb.mixt = 1, Niter = 10, discrete = 1, valid = 0, level = 0.05,
  newwindow = FALSE)

Arguments

times	vector of observation times
X	matrix of the M trajectories (each row is a trajectory with as much columns as observations)
model	name of the SDE: 'OU' (Ornstein-Uhlenbeck) or 'CIR' (Cox-Ingersoll-Ross)
drift.random	random effects in the drift: 0 if only fixed effects, 1 if one additive random effect, 2 if one multiplicative random effect or c(1,2) if 2 random effects. Default to c(1,2)
drift.fixed	NULL if the fixed effect(s) in the drift is (are) estimated, value of the fixed effect(s) otherwise. Default to NULL
diffusion.random	1 if σ is random, 0 otherwise. Default to 0
diffusion.fixed	NULL if σ is estimated (if fixed), value of σ otherwise. Default to NULL
nb.mixt	number of mixture components for the distribution of the random effects in the drift. Default to 1 (no mixture)
Niter	number of iterations for the EM algorithm if the random effects in the drift follow a mixture distribution. Default to 10
discrete	1 for using a contrast based on discrete observations, 0 otherwise. Default to 1
valid	1 if test validation, 0 otherwise. Default to 0
level	alpha for the predicion intervals. Default 0.05
newwindow	logical(1), if TRUE, a new window is opened for the plot. Default to FALSE

Details

Parametric estimation of the random effects density from M independent trajectories of the SDE:

dX_j(t)= (α_j- β_j X_j(t))dt + σ_j \ a(X_j(t)) dW_j(t),

j=1,…,M, where the (W_j(t)) are independant Wiener processes and the (X_j(t)) are observed without noise.

Specification of the random effects:

The drift includes no, one or two random effects:

1. if drift.random = 0: α_j \equiv α and β_j \equiv β are fixed

2. if drift.random = 1: β_j \equiv β is fixed and α_j is random

3. if drift.random = 2: α_j \equiv α is fixed and β_j is random

4. if drift.random = c(1,2): α_j and β_j are random

The diffusion includes either a fixed effect or a random effect:

1. if diffusion.random = 0: σ_j \equiv σ is fixed

2. if diffusion.random = 1: σ_j is random

Distribution of the random effects

If there is no random effect in the diffusion (diffusion.random = 0), there is at least on random effect in the drift that follows

1. a Gaussian distribution (nb.mixt=1): α_j \sim N(μ,Ω) or β_j \sim N(μ,Ω) or (α_j,β_j) \sim N(μ,Ω),

2. or a mixture of Gaussian distributions (nb.mixt=K, K>1): α_j \sim ∑_{k=1}^{K} p_k N(μ_k,Ω_k) or β_j \sim ∑_{k=1}^{K} p_k N(μ_k,Ω_k) or (α_j,β_j) \sim ∑_{k=1}^{K} p_k N(μ_k,Ω_k), where ∑_{k=1}^{K} p_k=1.

If there is one random effect in the diffusion (diffusion.random = 1), 1/σ_j^2 \sim Γ(a,λ), and the coefficients in the drift are conditionally Gaussian: α_j|σ_j \sim N(μ,σ_j^2 Ω) or β_j|σ_j \sim N(μ,σ_j^2 Ω) or (α_j,β_j)|σ_j \sim N(μ,σ_j^2 Ω), or they are fixed α_j \equiv α, β_j \equiv β.

SDEs

Two diffusions are implemented:

1. the Ornstein-Uhlenbeck model (OU) a(X_j(t))=1

2. the Cox-Ingersoll-Ross model (CIR) a(X_j(t))=√{X_j(t)}

Estimation

• If discrete = 0, the estimation is based on the exact likelihood associated with continuous observations ([1],[3]). This is only possible if diffusion.random = 0 and σ is not estimated by maximum likelihood but empirically by means of the quadratic variations.

• If discrete = 1, the likelihood of the Euler scheme of the mixed SDE is computed and maximized for estimating all the parameters.

• If nb.mixt > 1, an EM algorithm is implemented and the number of iterations of the algorithm must be specified with Niter.

• If valid = 1, two-thirds of the sample trajectories are used for estimation, while the rest is used for validation. A plot is then provided for visual comparison between the true trajectories of the test sample and some predicted trajectories simulated under the estimated model.

Value

index	is the vector of subscript in 1,...,M used for the estimation. Most of the time index=1:M, except for the CIR that requires positive trajectories.
estimphi	matrix of estimators of the drift random effects \hat{α}_j, or \hat{β}_j or (\hat{α}_j,\hat{β}_j)
estimpsi2	vector of estimators of the squared diffusion random effects \hat{σ}_j^2
gridf	grid of values for the plots of the random effects distribution in the drift, matrix form
gridg	grid of values for the plots of the random effects distribution in the diffusion, matrix form
estimf	estimator of the density of α_j, β_j or (α_j,β_j). Matrix form.
estimg	estimator of the density of σ_j^2. Matrix form.
mu	estimator of the mean of the random effects normal density
omega	estimator of the standard deviation of the random effects normal density
a	estimated value of the shape of the Gamma distribution
lambda	estimated value of the scale of the Gamma distribution
sigma2	value of the diffusion coefficient if it is fixed
bic	BIC criterium
aic	AIC criterium
model	initial choice
drift.random	initial choice
diffusion.random	initial choice
drift.fixed	initial choice
estim.drift.fix	1 if the fixed effects in the drift are estimated, 0 otherwise.
estim.diffusion.fixed	1 if the fixed effect in the diffusion is estimated, 0 otherwise.
discrete	initial choice
times	initial choice
X	initial choice

For mixture distributions in the drift:

mu	estimated value of the mean at each iteration of the algorithm. Niter x N x 2 array.
omega	estimated value of the standard deviation at each iteration of the algorithm. Niter x N x 2 array.
mixt.prop	estimated value of the mixture proportions at each iteration of the algorithm. Niter x N matrix.
probindi	posterior component probabilites. M x N matrix.

Author(s)

Maud Delattre and Charlotte Dion

References

See

[1] Maximum Likelihood Estimation for Stochastic Differential Equations with Random Effects, Delattre, M., Genon-Catalot, V. and Samson, A. Scandinavian Journal of Statistics 40(2) 2012 322-343

[2] Estimation of population parameters in stochastic differential equations with random effects in the diffusion coefficient, Delattre, M., Genon-Catalot, V. and Samson, A. ESAIM:PS 19 2015 671-688

[3] Mixtures of stochastic differential equations with random effects: application to data clustering, Delattre, M., Genon-Catalot, V. and Samson, A. Journal of Statistical Planning and Inference 173 2016 109-124

[4] Parametric inference for discrete observations of diffusion processes with mixed effects, Delattre, M., Genon-Catalot, V. and Laredo, C. hal-01332630 2016

[5] Estimation of the joint distribution of random effects for a discretely observed diffusion with random effects, Delattre, M., Genon-Catalot, V. and Laredo, C. hal-01446063 2017

Examples

# Example 1 : One random effect in the drift and one random effect in the diffusion

sim <- msde.sim(M = 25, T = 1, N = 1000, model = 'OU', 
                drift.random = 2, drift.param = c(0,0.5,0.5), 
                diffusion.random = 1, diffusion.param = c(8,1/2))

res <- msde.fit(times = sim$times, X = sim$X, model = 'OU', drift.random = 2, 
diffusion.random = 1)

summary(res)
plot(res)

## Not run: 

# Example 2 : one mixture of two random effects in the drift, and one fixed effect in
# the diffusion coefficient

sim <- msde.sim(M = 100, T = 5, N = 5000, model = 'OU', drift.random = c(1,2),
                diffusion.random = 0, 
                drift.param = matrix(c(0.5,1.8,0.25,0.25,1,2,0.25,0.25),nrow=2,byrow=FALSE), 
                diffusion.param = 0.1, nb.mixt = 2, mixt.prop = c(0.5,0.5))

# -- Estimation without validation
res <- msde.fit(times = sim$times, X = sim$X, model = 'OU', drift.random = c(1,2),
                nb.mixt=2, Niter = 10)

summary(res)
plot(res)

# -- Estimation with prediction
res.valid <- msde.fit(times = sim$times, X = sim$X, model = 'OU', drift.random = c(1,2),
                      nb.mixt=2, Niter = 10, valid = 1)

summary(res.valid)
plot(res.valid)

# Example 3 : CIR with one random effect in the drift and one random effect in the diffusion 
# coefficient

sim <- msde.sim(M = 100, T = 5, N = 5000, model = 'CIR', drift.random = 2,
                diffusion.random = 1, drift.param = c(4,1,0.1), diffusion.param = c(8,0.5),
                X0 = 1)

res <- msde.fit(times = sim$times, X = sim$X, model = 'CIR', drift.random = 2,
                diffusion.random = 1)

summary(res)

# Further examples can be found in the section dedicated to neuronal.data.
  
## End(Not run)

MsdeParEst documentation built on May 2, 2019, 3:01 a.m.