mixedsde.fit: Estimation Of The Random Effects In Mixed Stochastic...
In charlottedion/mixedsde: Estimation Methods for Stochastic Differential Mixed Effects Models
Description
Usage
Arguments
Details
Value
References
Examples
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
1
2
3
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, β 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.
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 φ=α, or β, or (α,β) 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 φ=α, or β 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 φ 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 φ, 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 φ, 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 α
beta
posterior samples (Markov chain) of β
mu
posterior samples (Markov chain) of μ
omega
posterior samples (Markov chain) of Ω
sigma2
posterior samples (Markov chain) of σ^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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
# 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)
charlottedion/mixedsde documentation built on May 13, 2019, 3:35 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Value References Examples
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
1 2 3 | 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, β 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.
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 φ=α, or β, or (α,β) from the efficient statitics (see |
estim.fixed |
if estim.fix is TRUE and random = 1 or 2, estimator of φ=α, or β from the efficient statitics (see |
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 φ 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 φ, 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 φ, 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 α |
beta |
posterior samples (Markov chain) of β |
mu |
posterior samples (Markov chain) of μ |
omega |
posterior samples (Markov chain) of Ω |
sigma2 |
posterior samples (Markov chain) of σ^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
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 | # 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)
|
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.