# 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
# 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)
Run the code above in your browser using DataLab