Estimation Of The Random Effects In Mixed Stochastic Differential Equations

Description

Estimation of the random effects (α_j, β_j) and of their density, parametrically or nonparametrically in the mixed SDE dX_j(t)= (α_j- β_j X_j(t))dt + σ a(X_j(t)) dW_j(t).

Usage

mixedsde.fit(times, X, model = c("OU", "CIR"), random, fixed = 0,
  estim.fix = 0, estim.method = c("nonparam", "paramML", "paramBayes"),
  gridf = NULL, prior, nMCMC = NULL)

Arguments

Details

Estimation of the random effects density from M independent trajectories of the SDE (the Brownian motions W_j are independent), with linear drift. Two diffusions are implemented, with one or two random effects:

Ornstein-Uhlenbeck model (OU)

If random = 1, β is a fixed effect: dX_j(t)= (α_j- β X_j(t))dt + σ dW_j(t)

If random = 2, α is a fixed effect: dX_j(t)= (α - β_j X_j(t))dt + σ dW_j(t)

If random = c(1,2), dX_j(t)= (α_j- β_j X_j(t))dt + σ dW_j(t)

Cox-Ingersoll-Ross model (CIR)

If random = 1, β is a fixed effect: dX_j(t)= (α_j- β X_j(t))dt + σ √{X_(t)} dWj_(t)

If random = 2, α is a fixed effect: dX_j(t)= (α - β_j X_j(t))dt + σ √{X_j(t)} dW_j(t)

If random = c(1,2), dX_j(t)= (α_j- β_j X_j(t))dt + σ √{X_j(t)} dW_j(t)

The nonparametric method estimates the density of the random effects with a kernel estimator (one-dimensional or two-dimensional density). The parametric method estimates the mean and standard deviation of the Gaussian distribution of the random effects.

Value

If there is a truncation threshold in the estimator

For the parametric maximum likelihood estimation

For the parametric Bayesian estimation

References

For the parametric estimation see: Maximum likelihood estimation for stochastic differential equations with random effects, M. Delattre, V. Genon-Catalot and A. Samson, Scandinavian Journal of Statistics 2012, Vol 40, 322–343

Bayesian Prediction of Crack Growth Based on a Hierarchical Diffusion Model. S. Hermann, K. Ickstadt and C. Mueller, appearing in: Applied Stochastic Models in Business and Industry 2016.

For the nonparametric estimation see:

Nonparametric estimation for stochastic differential equations with random effects, F. Comte, V. Genon-Catalot and A. Samson, Stochastic Processes and Their Applications 2013, Vol 7, 2522–2551

Estimation for stochastic differential equations with mixed effects, V. Genon-Catalot and C. Laredo 2014 e-print: hal-00807258

Bidimensional random effect estimation in mixed stochastic differential model, C. Dion and V. Genon-Catalot, Stochastic Inference for Stochastic Processes 2015, Springer Netherlands, 1–28

Examples

# Frequentist estimation
# Two random effects
model = 'CIR'; T <- 10
delta <- 0.1; M <- 100 # delta <- 0.001 and M <- 200 would yield good results
N <- floor(T/delta); sigma <- 0.01 ;
random <- c(1,2); density.phi <- 'gammainvgamma2'; param<- c(1.8, 0.8, 8, 0.05);
simu <- mixedsde.sim(M=M, T=T, N=N, model=model,random=random, density.phi=density.phi,
               param=param, sigma=sigma, invariant = 1)
X <- simu$X ; phi <- simu$phi; times <- simu$times
estim.method<- 'nonparam'
estim <- mixedsde.fit(times=times, X=X, model=model, random=random, estim.method= 'nonparam')
#To stock the results of the function, use method code{out}
#which put the outputs of the function on a list according to the frequentist or
# Bayesian approach.
outputsNP <- out(estim)

## Not run: 
plot(estim)
## End(Not run)
# It represents the bidimensional density, the histogram of the first estimated random
# effect eqn{alpha} with the  marginal of eqn{hat{f}} from the first coordonate which
# estimates  the density of eqn{alpha}. And the same for the second random effect
# eqn{eta}. More, it plots a qq-plot for the sample of estimator of the random effects
# compared with the quantiles of a Gaussian sample with the same mean and standard deviation.

summary(estim)
print(estim)
# Validation
# If numj is fixed by the user: this function simulates Mrep =100 (by default) new
# trajectories with the value of the estimated random effect. Then it plots on the
# left graph the Mrep new trajectories eqn{(Xnumj^{k}(t1), ... Xnumj^{k}(tN)),
# k= 1, ... Mrep} with in red the true trajectory eqn{(Xnumj(t1), ... Xnumj(tN))}.
#The right graph is a qq-plot of the quantiles of samples
# eqn{(Xnumj^{1}(ti), ... Xnumj^{Mrep}(ti))}
# for each time eqn{ti} compared with the uniform quantiles. The outputs of the function
# are: a matrix code{Xnew} dimension Mrepx N+1, vector of quantiles code{quantiles} length
# N and the number of the trajectory for the plot code{plotnumj= numj}
# If numj is not precised by the user, then, this function simulates Mrep =100 (by default)
# new trajectories for each estimated random effect. Then left graph is a plot of the Mrep
# new trajectories eqn{(Xj^{k}(t1), ... Xj^{k}(tN)), k= 1, ... Mrep}
#for a randomly chosen number j with in red the true trajectory eqn{(Xj(t1), ... Xj(tN))}.
#The right graph is a qq-plot of the quantiles of samples eqn{(Xj^{1}(ti), ... Xj^{Mrep}(ti))},
# for the same j and for each time eqn{ti}. The outputs of the function are: a list of
# matrices code{Xnew} length M, matrix of quantiles code{quantiles} dimension MxN
# and the number of the trajectory for the plot code{plotnumj}

validation <- valid(estim,  numj=floor(runif(1,1,M)))

# Parametric estimation
estim.method<-'paramML'
estim_param <- mixedsde.fit(times= times, X= X, model= model, random= random,
           estim.method = 'paramML')
outputsP <- out(estim_param)

#plot(estim_param)
summary(estim_param)

# Prediction for the frequentist approach
# This function uses the estimation of the density function to simulate a
# new sample of random effects according to this density. If code{plot.pred =1} (default)
# is plots on the top the predictive random effects versus the estimated random effects
# from the data. On the bottom, the left graph is the true trajectories, on the right
#the predictive trajectories and the empiric prediciton intervals at level
# code{level=0.05} (defaut). The function return on a list the prediction of phi
# code{phipred}, the prediction of X code{Xpred}, and the indexes of the
# corresponding true trajectories code{indexpred}

# Not run
## Not run: 
test1 <- pred(estim,  invariant  = 1)
test2 <- pred(estim_param, invariant  = 1)

## End(Not run)
# More graph
fhat <- outputsNP$estimf
fhat_trunc <- outputsNP$estimf.trunc
fhat_param <- outputsP$estimf

gridf <- outputsNP$gridf; gridf1 <- gridf[1,]; gridf2 <- gridf[2,]

marg1 <- ((max(gridf2)-min(gridf2))/length(gridf2))*apply(fhat,1,sum)
marg1_trunc <- ((max(gridf2)-min(gridf2))/length(gridf2))*apply(fhat_trunc,1,sum)
marg2 <- ((max(gridf1)-min(gridf1))/length(gridf1))*apply(fhat,2,sum)
marg2_trunc <- ((max(gridf1)-min(gridf1))/length(gridf1))*apply(fhat_trunc,2,sum)

marg1_param <- ((max(gridf2)-min(gridf2))/length(gridf2))*apply(fhat_param,1,sum)
marg2_param <- ((max(gridf1)-min(gridf1))/length(gridf1))*apply(fhat_param,2,sum)
f1 <-  (gridf1^(param[1]-1))*exp(-gridf1/param[2])/((param[2])^param[1]*gamma(param[1]))
f2 <-  (gridf2^(-param[3]-1)) * exp(-(1/param[4])*(1/gridf2)) *
 ((1/param[4])^param[3])*(1/gamma(param[3]))
par(mfrow=c(1,2))
plot(gridf1,f1,type='l', lwd=1,  xlab='', ylab='')
lines(gridf1,marg1_trunc,col='blue', lwd=2)
lines(gridf1,marg1,col='blue', lwd=2, lty=2)
lines(gridf1,marg1_param,col='red', lwd=2, lty=2)
plot(gridf2,f2,type='l', lwd=1, xlab='', ylab='')
lines(gridf2,marg2_trunc,col='blue', lwd=2)
lines(gridf2,marg2,col='blue', lwd=2, lty=2)
lines(gridf2,marg2_param,col='red', lwd=2, lty=2)

cutoff <- outputsNP$cutoff
phihat <- outputsNP$estimphi
phihat_trunc <- outputsNP$estimphi.trunc
par(mfrow=c(1,2))
plot.ts(phi[1,], phihat[1,], xlim=c(0, 15), ylim=c(0,15), pch=18); abline(0,1)
points(phi[1,]*(1-cutoff), phihat[1,]*(1-cutoff), xlim=c(0, 20), ylim=c(0,20),pch=18, col='red');
abline(0,1)
plot.ts(phi[2,], phihat[2,], xlim=c(0, 15), ylim=c(0,15),pch=18); abline(0,1)
points(phi[2,]*(1-cutoff), phihat[2,]*(1-cutoff), xlim=c(0, 20), ylim=c(0,20),pch=18, col='red');
abline(0,1)

# one random effect:
## Not run: 
model <-'OU'
random <- 1
M <- 80; T <- 100  ; N <- 2000
sigma <- 0.1 ; X0 <- 0
density.phi <- 'normal'
fixed <- 2 ; param <- c(1, 0.2)
#-------------------
#- simulation
simu <- mixedsde.sim(M,T=T,N=N,model=model,random=random, fixed=fixed,density.phi=density.phi,
               param=param, sigma=sigma, X0=X0)
X <- simu$X
phi <- simu$phi
times <-simu$times
plot(times, X[10,], type='l')

#- parametric estimation
estim.method<-'paramML'
estim_param <- mixedsde.fit(times, X=X, model=model, random=random, estim.fix= 1,
               estim.method=estim.method)
outputsP <- out(estim_param)
estim.fixed <- outputsP$estim.fixed #to compare with fixed
#or
estim_param2 <- mixedsde.fit(times, X=X, model=model, random=random, fixed = fixed,
             estim.method=estim.method)
outputsP2 <- out(estim_param2)
#- nonparametric estimation
estim.method <- 'nonparam'
estim <- mixedsde.fit(times, X, model=model, random=random, fixed = fixed,
           estim.method=estim.method)
outputsNP <- out(estim)

plot(estim)
print(estim)
summary(estim)

plot(estim_param)
print(estim_param)
summary(estim_param)

valid1 <- valid(estim,  numj=floor(runif(1,1,M)))
test1 <- pred(estim )
test2 <- pred(estim_param)

## End(Not run)

# Parametric Bayesian estimation
# one random effect
random <- 1; sigma <- 0.1; fixed <- 5; param <- c(3, 0.5)
sim <- mixedsde.sim(M = 20, T = 1, N = 50, model = 'OU', random = random, fixed = fixed,
       density.phi = 'normal',param= param, sigma= sigma, X0 = 0, op.plot = 1)

# here: only 100 iterations for example - should be much more!
prior <- list( m = c(param[1], fixed), v = c(param[1], fixed), alpha.omega = 11,
            beta.omega = param[2]^2*10, alpha.sigma = 10, beta.sigma = sigma^2*9)
estim_Bayes <- mixedsde.fit(times = sim$times, X = sim$X, model = 'OU', random,
           estim.method = 'paramBayes', prior = prior, nMCMC = 100)

validation <- valid(estim_Bayes, numj = 10)
plot(estim_Bayes)
outputBayes <- out(estim_Bayes)
summary(outputBayes)
(results_Bayes <- summary(estim_Bayes))
plot(estim_Bayes, style = 'cred.int', true.phi = sim$phi)
print(estim_Bayes)
## Not run: 
pred.result <- pred(estim_Bayes)
summary(out(pred.result))
plot(pred.result)

pred.result.trajectories <- pred(estim_Bayes, trajectories = TRUE)

## End(Not run)
# second example
## Not run: 
random <- 2; sigma <- 0.2; fixed <- 5; param <- c(3, 0.5)
sim <- mixedsde.sim(M = 20, T = 1, N = 100, model = 'CIR', random = random,
        fixed = fixed, density.phi = 'normal',param = param, sigma = sigma, X0 = 0.1, op.plot = 1)

prior <- list(m = c(fixed, param[1]), v = c(fixed, param[1]), alpha.omega = 11,
         beta.omega = param[2]^2*10, alpha.sigma = 10, beta.sigma = sigma^2*9)

estim_Bayes <- mixedsde.fit(times = sim$times, X = sim$X, model = 'CIR', random = random,
                 estim.method = 'paramBayes', prior = prior, nMCMC = 1000)

pred.result <- pred(estim_Bayes)

## End(Not run)



# invariant case
## Not run: 
random <- 1; sigma <- 0.1; fixed <- 5; param <- c(3, 0.5)
sim <- mixedsde.sim(M = 50, T = 5, N = 100, model = 'OU', random = random, fixed = fixed,
           density.phi = 'normal',param = param, sigma = sigma, invariant = 1, op.plot = 1)

prior <- list(m = c(param[1], fixed), v = c(param[1], 1e-05), alpha.omega = 11,
       beta.omega = param[2]^2*10, alpha.sigma = 10, beta.sigma = sigma^2*9)
estim_Bayes <- mixedsde.fit(times = sim$times, X = sim$X, model = 'OU', random,
       estim.method = 'paramBayes', prior = prior, nMCMC = 100)
plot(estim_Bayes)

pred.result <- pred(estim_Bayes, invariant = 1)
pred.result.traj <- pred(estim_Bayes, invariant = 1, trajectories = TRUE)

## End(Not run)

Results