mixedsde.fit {mixedsde}

Estimation Of The Random Effects In Mixed Stochastic Differential Equations

Description

Estimation of the random effects (\alpha_j, \beta_j) and of their density, parametrically or nonparametrically in the mixed SDE dX_j(t)= (\alpha_j- \beta_j X_j(t))dt + \sigma 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

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)

random

random effects in the drift: 1 if one additive random effect, 2 if one multiplicative random effect or c(1,2) if 2 random effects

fixed

fixed effect in the drift: value of the fixed effect when there is only one random effect and it is not estimated, 0 otherwise

estim.fix

default 0, 1 if random = 1 or 2, method = 'paramML' then the fixed parameter is estimated

estim.method

estimation method: 'paramML' for a Gaussian parametric estimation by maximum likelihood, 'paramBayes' for a Gaussian parametric Bayesian estimation or 'nonparam' for a non-parametric estimation

gridf

if nonparametric estimation: grid of values on which the density is estimated, a matrix with two rows if two random effects; NULL by default and then grid is chosen as a function of the estimated values of the random effects. For the plots this grid is used.

prior

if method = 'paramBayes', list of prior parameters: mean and variance of the Gaussian prior on the mean mu, shape and scale of the inverse Gamma prior for the variances omega, shape and scale of the inverse Gamma prior for sigma

nMCMC

if method = 'paramBayes', number of iterations of the MCMC algorithm

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, \beta is a fixed effect: dX_j(t)= (\alpha_j- \beta X_j(t))dt + \sigma dW_j(t)

If random = 2, \alpha is a fixed effect: dX_j(t)= (\alpha - \beta_j X_j(t))dt + \sigma dW_j(t)

If random = c(1,2), dX_j(t)= (\alpha_j- \beta_j X_j(t))dt + \sigma dW_j(t)

Cox-Ingersoll-Ross model (CIR)

If random = 1, \beta is a fixed effect: dX_j(t)= (\alpha_j- \beta X_j(t))dt + \sigma \sqrt{X_(t)} dWj_(t)

If random = 2, \alpha is a fixed effect: dX_j(t)= (\alpha - \beta_j X_j(t))dt + \sigma \sqrt{X_j(t)} dW_j(t)

If random = c(1,2), dX_j(t)= (\alpha_j- \beta_j X_j(t))dt + \sigma \sqrt{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.

Validation method: For a number of trajectory numj (fixed by the user or randomly chosen) 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 (Xnumj^{k}(t1), ... Xnumj^{k}(tN)), k= 1, ... Mrep with in red the true trajectory (Xnumj(t1), ... Xnumj(tN)). The right graph is a qq-plot of the quantiles of samples (Xnumj^{1}(ti), ... Xnumj^{Mrep}(ti)) for each time ti compared with the uniform quantiles. The outputs of the function are: a matrix Xnew dimension Mrepx N+1, vector of quantiles quantiles length N and the number of the trajectory for the plot numj

Prediction method 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 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 level=0.05 (defaut). The function return on a list the prediction of phi phipred, the prediction of X Xpred, and the indexes of the corresponding true trajectories indexpred

Value

index

is the vector of subscript in 1,...,M where the estimation of phi has been done, most of the time index= 1:M

estimphi

matrix of estimators of \phi=\alpha, or \beta, or (\alpha,\beta) from the efficient statitics (see UV), matrix of two lines if random =c(1,2), numerical type otherwise

estim.fixed

if estim.fix is TRUE and random = 1 or 2, estimator of \phi=\alpha, or \beta from the efficient statitics (see UV), 0 otherwise

gridf

grid of values on which the estimated is done for the nonparametric method, otherwise, grid used for the plots, matrix form

estimf

estimator of the density of \phi from a kernel estimator from package: stats, function: density, or package: MASS, function: kde2D. Matrix form: one line if one random effect or square matrix otherwise

If there is a truncation threshold in the estimator

cutoff

the binary vector of cutoff, FALSE otherwise

estimphi.trunc

troncated estimator of \phi, vector or matrix of 0 if we do not use truncation, matrix of two lines if random =c(1,2), numerical type otherwise

estimf.trunc

troncated estimator of the density of \phi, vector or matrix of 0 if we do not use truncation, matrix if random =c(1,2), numerical type otherwise

For the parametric maximum likelihood estimation

mu

estimator of the mean of the random effects normal density, 0 if we do nonparametric estimation

omega

estimator of the standard deviation of the random effects normal density, 0 if we do nonparametric estimation

bic

BIC criterium, 0 if we do nonparametric estimation

aic

AIC criterium, 0 if we do nonparametric estimation

model

initial choice

random

initial choice

fixed

initial choice

times

initial choice

X

initial choice

For the parametric Bayesian estimation

alpha

posterior samples (Markov chain) of \alpha

beta

posterior samples (Markov chain) of \beta

mu

posterior samples (Markov chain) of \mu

omega

posterior samples (Markov chain) of \Omega

sigma2

posterior samples (Markov chain) of \sigma^2

model

initial choice

random

initial choice

burnIn

proposal for burn-in period

thinning

proposal for thinning rate

prior

initial choice or calculated by the first 10% of series

times

initial choice

X

initial choice

ind.4.prior

in the case of calculation of prior parameters: the indices of used series

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{\beta}. 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 

validation <- valid(estim)

# 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)

# 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)
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)