Package 'yuima'

Description

The adabayes.mcmc class is a class of the yuima package that extends the mle-class.

Usage

adaBayes(yuima, start, prior, lower, upper, fixed = numeric(0), envir = globalenv(),
method = "mcmc", iteration = NULL,mcmc,rate =1, rcpp = TRUE, 
algorithm = "randomwalk",center=NULL,sd=NULL,rho=NULL,
path = FALSE)

Arguments

Details

Calculate the Bayes estimator for stochastic processes by using the quasi-likelihood function. The calculation is performed by the Markov chain Monte Carlo method. Currently, the Random-walk Metropolis algorithm and the Mixed preconditioned Crank-Nicolson algorithm is implemented.

Slots

mcmc:

is a list of MCMC objects for all estimated parameters.

accept_rate:

is a list acceptance rates for diffusion and drift parts.

call:

is an object of class language.

fullcoef:

is an object of class list that contains estimated parameters.

vcov:

is an object of class matrix.

coefficients:

is an object of class vector that contains estimated parameters.

Note

algorithm = nomcmc is unstable.

Author(s)

Kengo Kamatani with YUIMA project Team

References

Yoshida, N. (2011). Polynomial type large deviation inequalities and quasi-likelihood analysis for stochastic differential equations. Annals of the Institute of Statistical Mathematics, 63(3), 431-479. Uchida, M., & Yoshida, N. (2014). Adaptive Bayes type estimators of ergodic diffusion processes from discrete observations. Statistical Inference for Stochastic Processes, 17(2), 181-219. Kamatani, K. (2017). Ergodicity of Markov chain Monte Carlo with reversible proposal. Journal of Applied Probability, 54(2).

Examples

## Not run: 
set.seed(123)
b <- c("-theta1*x1+theta2*sin(x2)+50","-theta3*x2+theta4*cos(x1)+25")
a <- matrix(c("4+theta5","1","1","2+theta6"),2,2)
true = list(theta1 = 0.5, theta2 = 5,theta3 = 0.3,
            theta4 = 5, theta5 = 1, theta6 = 1)
lower = list(theta1=0.1,theta2=0.1,theta3=0,
             theta4=0.1,theta5=0.1,theta6=0.1)
upper = list(theta1=1,theta2=10,theta3=0.9,
             theta4=10,theta5=10,theta6=10)
start = list(theta1=runif(1),
             theta2=rnorm(1),
             theta3=rbeta(1,1,1),
             theta4=rnorm(1),
             theta5=rgamma(1,1,1),
             theta6=rexp(1))
yuimamodel <- setModel(drift=b,diffusion=a,state.variable=c("x1", "x2"),solve.variable=c("x1","x2"))
yuimasamp <- setSampling(Terminal=50,n=50*10)
yuima <- setYuima(model = yuimamodel, sampling = yuimasamp)
yuima <- simulate(yuima, xinit = c(100,80),
                  true.parameter = true,sampling = yuimasamp)
prior <-
  list(
    theta1=list(measure.type="code",df="dunif(z,0,1)"),
    theta2=list(measure.type="code",df="dnorm(z,0,1)"),
    theta3=list(measure.type="code",df="dbeta(z,1,1)"),
    theta4=list(measure.type="code",df="dgamma(z,1,1)"),
    theta5=list(measure.type="code",df="dnorm(z,0,1)"),
    theta6=list(measure.type="code",df="dnorm(z,0,1)")
  )
set.seed(123)
mle <- qmle(yuima, start = start, lower = lower, upper = upper, method = "L-BFGS-B",rcpp=TRUE)
print(mle@coef)
center<-list(theta1=0.5,theta2=5,theta3=0.3,theta4=4,theta5=3,theta6=3)
sd<-list(theta1=0.001,theta2=0.001,theta3=0.001,theta4=0.01,theta5=0.5,theta6=0.5)
bayes <- adaBayes(yuima, start=start, prior=prior,lower=lower,upper=upper,
                  method="mcmc",mcmc=1000,rate = 1, rcpp = TRUE,
                   algorithm = "randomwalk",center = center,sd=sd,
                   path=TRUE)
print(bayes@fullcoef)
print(bayes@accept_rate)
print(bayes@mcmc$theta1[1:10])

## End(Not run)

Description

Asymptotic expansion of uni-dimensional and multi-dimensional diffusion processes.

Usage

ae(
  model,
  xinit,
  order = 1L,
  true.parameter = list(),
  sampling = NULL,
  eps.var = "eps",
  solver = "rk4",
  verbose = FALSE
)

Arguments

Details

If sampling is not provided, then model must be an object of yuima-class with non-empty sampling.

if eps.var does not appear in the model specification, then it is internally added in front of the diffusion matrix to apply the asymptotic expansion scheme.

Value

An object of yuima.ae-class

Author(s)

Emanuele Guidotti <[email protected]>

Examples

## Not run: 
# model
gbm <- setModel(drift = 'mu*x', diffusion = 'sigma*x', solve.variable = 'x')

# settings
xinit <- 100
par <- list(mu = 0.01, sigma = 0.2)
sampling <- setSampling(Initial = 0, Terminal = 1, n = 1000)

# asymptotic expansion
approx <- ae(model = gbm, sampling = sampling, order = 4, true.parameter = par, xinit = xinit)

# exact density
x <- seq(50, 200, by = 0.1)
exact <- dlnorm(x = x, meanlog = log(xinit)+(par$mu-0.5*par$sigma^2)*1, sdlog = par$sigma*sqrt(1))

# compare
plot(x, exact, type = 'l', ylab = "Density")
lines(x, aeDensity(x = x, ae = approx, order = 1), col = 2)
lines(x, aeDensity(x = x, ae = approx, order = 2), col = 3)
lines(x, aeDensity(x = x, ae = approx, order = 3), col = 4)
lines(x, aeDensity(x = x, ae = approx, order = 4), col = 5)

## End(Not run)

Description

Asymptotic Expansion - Characteristic Function

Usage

aeCharacteristic(..., ae, eps = 1, order = NULL)

Arguments

Value

Characteristic function evaluated on the given grid.

Examples

## Not run: 
# model
gbm <- setModel(drift = 'mu*x', diffusion = 'sigma*x', solve.variable = 'x')

# settings
xinit <- 100
par <- list(mu = 0.01, sigma = 0.2)
sampling <- setSampling(Initial = 0, Terminal = 1, n = 1000)

# asymptotic expansion
approx <- ae(model = gbm, sampling = sampling, order = 4, true.parameter = par, xinit = xinit)

# The following are all equivalent methods to specify the grid via ....
# Notice that the character 'u1' corresponds to the 'u.var' of the ae object.
approx@u.var

# 1) named argument  
u1 <- seq(0, 1, by = 0.1)
psi <- aeCharacteristic(u1 = u1, ae = approx, order = 4)
# 2) data frame
df <- data.frame(u1 = seq(0, 1, by = 0.1))
psi <- aeCharacteristic(df, ae = approx, order = 4)
# 3) environment
env <- new.env()
env$u1 <- seq(0, 1, by = 0.1)
psi <- aeCharacteristic(env, ae = approx, order = 4)
# 4) list
lst <- list(u1 = seq(0, 1, by = 0.1))
psi <- aeCharacteristic(lst, ae = approx, order = 4)

## End(Not run)

Description

Asymptotic Expansion - Density

Usage

aeDensity(..., ae, eps = 1, order = NULL)

Arguments

Value

Probability density function evaluated on the given grid.

Examples

## Not run: 
# model
gbm <- setModel(drift = 'mu*x', diffusion = 'sigma*x', solve.variable = 'x')

# settings
xinit <- 100
par <- list(mu = 0.01, sigma = 0.2)
sampling <- setSampling(Initial = 0, Terminal = 1, n = 1000)

# asymptotic expansion
approx <- ae(model = gbm, sampling = sampling, order = 4, true.parameter = par, xinit = xinit)

# The following are all equivalent methods to specify the grid via ....
# Notice that the character 'x' corresponds to the solve.variable of the yuima model.

# 1) named argument  
x <- seq(50, 200, by = 0.1)
density <- aeDensity(x = x, ae = approx, order = 4)
# 2) data frame
df <- data.frame(x = seq(50, 200, by = 0.1))
density <- aeDensity(df, ae = approx, order = 4)
# 3) environment
env <- new.env()
env$x <- seq(50, 200, by = 0.1)
density <- aeDensity(env, ae = approx, order = 4)
# 4) list
lst <- list(x = seq(50, 200, by = 0.1))
density <- aeDensity(lst, ae = approx, order = 4)

# exact density
exact <- dlnorm(x = x, meanlog = log(xinit)+(par$mu-0.5*par$sigma^2)*1, sdlog = par$sigma*sqrt(1))

# compare
plot(x = exact, y = density, xlab = "Exact", ylab = "Approximated")

## End(Not run)

Description

Compute the expected value of functionals.

Usage

aeExpectation(f, bounds, ae, eps = 1, order = NULL, ...)

Arguments

Value

return value of cubintegrate. The expectation of the functional provided.

Examples

## Not run: 
# model
gbm <- setModel(drift = 'mu*x', diffusion = 'sigma*x', solve.variable = 'x')

# settings
xinit <- 100
par <- list(mu = 0.01, sigma = 0.2)
sampling <- setSampling(Initial = 0, Terminal = 1, n = 1000)

# asymptotic expansion
approx <- ae(model = gbm, sampling = sampling, order = 4, true.parameter = par, xinit = xinit)

# compute the mean via integration
aeExpectation(f = 'x', bounds = list(x = c(0,1000)), ae = approx)

# compare with the mean computed by differentiation of the characteristic function
aeMean(approx)

## End(Not run)

Description

Asymptotic Expansion - Kurtosis

Usage

aeKurtosis(ae, eps = 1, order = NULL)

Arguments

Value

numeric.

Examples

## Not run: 
# model
gbm <- setModel(drift = 'mu*x', diffusion = 'sigma*x', solve.variable = 'x')

# settings
xinit <- 100
par <- list(mu = 0.01, sigma = 0.2)
sampling <- setSampling(Initial = 0, Terminal = 1, n = 1000)

# asymptotic expansion
approx <- ae(model = gbm, sampling = sampling, order = 4, true.parameter = par, xinit = xinit)

# expansion order max
aeKurtosis(ae = approx)

# expansion order 1
aeKurtosis(ae = approx, order = 1)

## End(Not run)

Description

Asymptotic Expansion - Marginals

Usage

aeMarginal(ae, var)

Arguments

Value

An object of yuima.ae-class

Examples

## Not run: 
# multidimensional model
gbm <- setModel(drift = c('mu*x1','mu*x2'), 
                diffusion = matrix(c('sigma1*x1',0,0,'sigma2*x2'), nrow = 2), 
                solve.variable = c('x1','x2'))

# settings
xinit <- c(100, 100)
par <- list(mu = 0.01, sigma1 = 0.2, sigma2 = 0.1)
sampling <- setSampling(Initial = 0, Terminal = 1, n = 1000)

# asymptotic expansion
approx <- ae(model = gbm, sampling = sampling, order = 3, true.parameter = par, xinit = xinit)

# extract marginals
margin1 <- aeMarginal(ae = approx, var = "x1")
margin2 <- aeMarginal(ae = approx, var = "x2")

# compare with exact solution for marginal 1
x1 <- seq(50, 200, by = 0.1)
exact <- dlnorm(x = x1, meanlog = log(xinit[1])+(par$mu-0.5*par$sigma1^2), sdlog = par$sigma1)
plot(x1, exact, type = 'p', ylab = "Density")
lines(x1, aeDensity(x1 = x1, ae = margin1, order = 3), col = 2)

# compare with exact solution for marginal 2
x2 <- seq(50, 200, by = 0.1)
exact <- dlnorm(x = x2, meanlog = log(xinit[2])+(par$mu-0.5*par$sigma2^2), sdlog = par$sigma2)
plot(x2, exact, type = 'p', ylab = "Density")
lines(x2, aeDensity(x2 = x2, ae = margin2, order = 3), col = 2)


## End(Not run)

Description

Asymptotic Expansion - Mean

Usage

aeMean(ae, eps = 1, order = NULL)

Arguments

Value

numeric.

Examples

## Not run: 
# model
gbm <- setModel(drift = 'mu*x', diffusion = 'sigma*x', solve.variable = 'x')

# settings
xinit <- 100
par <- list(mu = 0.01, sigma = 0.2)
sampling <- setSampling(Initial = 0, Terminal = 1, n = 1000)

# asymptotic expansion
approx <- ae(model = gbm, sampling = sampling, order = 4, true.parameter = par, xinit = xinit)

# expansion order max
aeMean(ae = approx)

# expansion order 1
aeMean(ae = approx, order = 1)

## End(Not run)

Description

Asymptotic Expansion - Moments

Usage

aeMoment(ae, m = 1, eps = 1, order = NULL)

Arguments

Value

numeric.

Examples

## Not run: 
# model
gbm <- setModel(drift = 'mu*x', diffusion = 'sigma*x', solve.variable = 'x')

# settings
xinit <- 100
par <- list(mu = 0.01, sigma = 0.2)
sampling <- setSampling(Initial = 0, Terminal = 1, n = 1000)

# asymptotic expansion
approx <- ae(model = gbm, sampling = sampling, order = 4, true.parameter = par, xinit = xinit)

# second moment, expansion order max
aeMoment(ae = approx, m = 2)

# second moment, expansion order 3
aeMoment(ae = approx, m = 2, order = 3)

# second moment, expansion order 2
aeMoment(ae = approx, m = 2, order = 2)

# second moment, expansion order 1
aeMoment(ae = approx, m = 2, order = 1)

## End(Not run)

Description

Asymptotic Expansion - Standard Deviation

Usage

aeSd(ae, eps = 1, order = NULL)

Arguments

Value

numeric.

Examples

## Not run: 
# model
gbm <- setModel(drift = 'mu*x', diffusion = 'sigma*x', solve.variable = 'x')

# settings
xinit <- 100
par <- list(mu = 0.01, sigma = 0.2)
sampling <- setSampling(Initial = 0, Terminal = 1, n = 1000)

# asymptotic expansion
approx <- ae(model = gbm, sampling = sampling, order = 4, true.parameter = par, xinit = xinit)

# expansion order max
aeSd(ae = approx)

# expansion order 1
aeSd(ae = approx, order = 1)

## End(Not run)

Description

Asymptotic Expansion - Skewness

Usage

aeSkewness(ae, eps = 1, order = NULL)

Arguments

Value

numeric.

Examples

## Not run: 
# model
gbm <- setModel(drift = 'mu*x', diffusion = 'sigma*x', solve.variable = 'x')

# settings
xinit <- 100
par <- list(mu = 0.01, sigma = 0.2)
sampling <- setSampling(Initial = 0, Terminal = 1, n = 1000)

# asymptotic expansion
approx <- ae(model = gbm, sampling = sampling, order = 4, true.parameter = par, xinit = xinit)

# expansion order max
aeSkewness(ae = approx)

# expansion order 1
aeSkewness(ae = approx, order = 1)

## End(Not run)

Description

calculate the fisrt and second term of asymptotic expansion of the functional mean.

Usage

asymptotic_term(yuima, block=100, rho, g, expand.var="e")

Arguments

Details

Calculate the first and second term of asymptotic expansion of the expected value of the functional associated with a sde. The returned value d0 + epsilon * d1 is approximation of the expected value.

Value

Note

we need to fix this routine.

Author(s)

YUIMA Project Team

Examples

## Not run: 
# to the Black-Scholes economy:
# dXt^e = Xt^e * dt + e * Xt^e * dWt
diff.matrix <- "x*e"
model <- setModel(drift = "x", diffusion = diff.matrix)
# call option is evaluated by averating
# max{ (1/T)*int_0^T Xt^e dt, 0}, the first argument is the functional of interest:
Terminal <- 1
xinit <- c(1)
f <- list( c(expression(x/Terminal)), c(expression(0)))
F <- 0
division <- 1000
e <- .3
yuima <- setYuima(model = model, sampling = setSampling(Terminal=Terminal, n=division))
yuima <- setFunctional( yuima, f=f,F=F, xinit=xinit,e=e)

# asymptotic expansion
rho <- expression(0)
F0 <- F0(yuima)
get_ge <- function(x,epsilon,K,F0){
  tmp <- (F0 - K) + (epsilon * x) 
  tmp[(epsilon * x) < (K-F0)] <- 0
  return( tmp )
}
g <- function(x) get_ge(x,epsilon=e,K=1,F0=F0)
set.seed(123)
asymp <- asymptotic_term(yuima, block=10, rho,g)
asymp
sum(asymp$d0 + e * asymp$d1)


### An example of multivariate case: Heston model
## a <- 1;C <- 1;d <- 10;R<-.1
## diff.matrix <- matrix( c("x1*sqrt(x2)*e", "e*R*sqrt(x2)",0,"sqrt(x2*(1-R^2))*e"), 2,2)
## model <- setModel(drift = c("a*x1","C*(10-x2)"), 
## diffusion = diff.matrix,solve.variable=c("x1","x2"),state.variable=c("x1","x2"))
## call option is evaluated by averating
## max{ (1/T)*int_0^T Xt^e dt, 0}, the first argument is the functional of interest:
##
## Terminal <- 1
## xinit <- c(1,1)
## 
## f <- list( c(expression(0), expression(0)),   
## c(expression(0), expression(0)) , c(expression(0), expression(0))  )
## F <- expression(x1,x2)
## 
## division <- 1000
## e <- .3
## 
## yuima <- setYuima(model = model, sampling = setSampling(Terminal=Terminal, n=division))
## yuima <- setFunctional( yuima, f=f,F=F, xinit=xinit,e=e)
## 
## rho <- expression(x1)
## F0 <- F0(yuima)
## get_ge <- function(x){
##  return( max(x[1],0))
## }
## g <- function(x) get_ge(x)
## set.seed(123)
## asymp <- asymptotic_term(yuima, block=10, rho,g)
## sum(asymp$d0 + e * asymp$d1)


## End(Not run)

Description

Tests the presence of jumps using the statistic proposed in Barndorff-Nielsen and Shephard (2004,2006) for each component.

Usage

bns.test(yuima, r = rep(1, 4), type = "standard", adj = TRUE)

Arguments

Details

For the i-th component, the test statistic is equal to the i-th component of sqrt(n)*(mpv(yuima,2)-mpv(yuima,c(1,1)))/sqrt(vartheta*mpv(yuima,r)) when type="standard", sqrt(n)*log(mpv(yuima,2)/mpv(yuima,c(1,1)))/sqrt(vartheta*mpv(yuima,r)/mpv(yuima,c(1,1))^2) when type="log" and sqrt(n)*(1-mpv(yuima,c(1,1))/mpv(yuima,2))/sqrt(vartheta*mpv(yuima,r)/mpv(yuima,c(1,1))^2) when type="ratio". Here, n is equal to the length of the i-th component of the zoo.data of yuima minus 1 and vartheta is pi^2/4+pi-5. When adj=TRUE, mpv(yuima,r)[i]/mpv(yuima,c(1,1))^2)[i] is replaced with 1 if it is less than 1.

Value

A list with the same length as the zoo.data of yuima. Each component of the list has class “htest” and contains the following components:

Note

Theoretically, this test may be invalid if sampling is irregular.

Author(s)

Yuta Koike with YUIMA Project Team

References

Barndorff-Nielsen, O. E. and Shephard, N. (2004) Power and bipower variation with stochastic volatility and jumps, Journal of Financial Econometrics, 2, no. 1, 1–37.

Barndorff-Nielsen, O. E. and Shephard, N. (2006) Econometrics of testing for jumps in financial economics using bipower variation, Journal of Financial Econometrics, 4, no. 1, 1–30.

Huang, X. and Tauchen, G. (2005) The relative contribution of jumps to total price variance, Journal of Financial Econometrics, 3, no. 4, 456–499.

See Also

lm.jumptest, mpv, minrv.test, medrv.test, pz.test

Examples

set.seed(123)

# One-dimensional case
## Model: dXt=t*dWt+t*dzt, 
## where zt is a compound Poisson process with intensity 5 and jump sizes distribution N(0,0.1).

model <- setModel(drift=0,diffusion="t",jump.coeff="t",measure.type="CP",
                  measure=list(intensity=5,df=list("dnorm(z,0,sqrt(0.1))")),
                  time.variable="t")

yuima.samp <- setSampling(Terminal = 1, n = 390) 
yuima <- setYuima(model = model, sampling = yuima.samp) 
yuima <- simulate(yuima)
plot(yuima) # The path seems to involve some jumps

bns.test(yuima) # standard type

bns.test(yuima,type="log") # log type

bns.test(yuima,type="ratio") # ratio type

# Multi-dimensional case
## Model: dXkt=t*dWk_t (k=1,2,3) (no jump case).

diff.matrix <- diag(3)
diag(diff.matrix) <- c("t","t","t")
model <- setModel(drift=c(0,0,0),diffusion=diff.matrix,time.variable="t",
                  solve.variable=c("x1","x2","x3"))

yuima.samp <- setSampling(Terminal = 1, n = 390) 
yuima <- setYuima(model = model, sampling = yuima.samp) 
yuima <- simulate(yuima)
plot(yuima)

bns.test(yuima)

Description

The carma.info-class is a class of the yuima package.

Details

The carma.info-class object cannot be directly specified by the user but it is constructed when the yuima.carma-class object is constructed via setCarma.

Slots

p:

Number of autoregressive coefficients.

q:

Number of moving average coefficients.

loc.par:

Label of location coefficient.

scale.par:

Label of scale coefficient.

ar.par:

Label of autoregressive coefficients.

ma.par:

Label of moving average coefficients.

lin.par:

Label of linear coefficients.

Carma.var:

Label of the observed process.

Latent.var:

Label of the unobserved process.

XinExpr:

Logical variable. If XinExpr=FALSE, the starting condition of Latent.var is zero otherwise each component of Latent.var has a parameter as a starting point.

Author(s)

The YUIMA Project Team

Description

The carmaHawkes.info-class is a class of the yuima package.

Details

The carmaHawkes.info-class object cannot be directly specified by the user but it is constructed when the yuima.carmaHawkes-class object is constructed via setCarmaHawkes.

Slots

p:

Number of autoregressive coefficients.

q:

Number of moving average coefficients.

Counting.Process:

Label of Counting process.

base.Int:

Label of baseline Intensity parameter.

ar.par:

Label of autoregressive coefficients.

ma.par:

Label of moving average coefficients.

Intensity.var:

Label of the Intensity process.

Latent.var:

Label of the unobserved process.

XinExpr:

Logical variable. If XinExpr=FALSE, the starting condition of Latent.var is zero otherwise each component of Latent.var has a parameter as a starting point.

Type.Jump:

Logical variable. If XinExpr=TRUE, the jump size is deterministic

Author(s)

The YUIMA Project Team

Contacts: Lorenzo Mercuri [email protected]

Description

Retrieve the increment of the underlying Levy for the carma(p,q) process using the approach developed in Brockwell et al.(2011)

Usage

CarmaNoise(yuima, param, data=NULL, NoNeg.Noise=FALSE)

Arguments

Value

Note

The function qmle uses the function CarmaNoise for estimation of underlying Levy in the carma model.

Author(s)

The YUIMA Project Team

References

Brockwell, P., Davis, A. R. and Yang. Y. (2011) Estimation for Non-Negative Levy-Driven CARMA Process, Journal of Business And Economic Statistics, 29 - 2, 250-259.

Examples

## Not run: 
#Ex.1: Carma(p=3, q=0) process driven by a brownian motion.

mod0<-setCarma(p=3,q=0)

# We fix the autoregressive and moving average parameters
# to ensure the existence of a second order stationary solution for the process.

true.parm0 <-list(a1=4,a2=4.75,a3=1.5,b0=1)

# We simulate a trajectory of the Carma model.

numb.sim<-1000
samp0<-setSampling(Terminal=100,n=numb.sim)
set.seed(100)
incr.W<-matrix(rnorm(n=numb.sim,mean=0,sd=sqrt(100/numb.sim)),1,numb.sim)

sim0<-simulate(mod0,
               true.parameter=true.parm0,
               sampling=samp0, increment.W=incr.W)

#Applying the CarmaNoise

system.time(
  inc.Levy0<-CarmaNoise(sim0,true.parm0)
)

# We compare the orginal with the estimated noise increments 

par(mfrow=c(1,2))
plot(t(incr.W)[1:998],type="l", ylab="",xlab="time")
title(main="True Brownian Motion",font.main="1")
plot(inc.Levy0,type="l", main="Filtered Brownian Motion",font.main="1",ylab="",xlab="time")

# Ex.2: carma(2,1) driven  by a compound poisson
# where jump size is normally distributed and
# the lambda is equal to 1.

mod1<-setCarma(p=2,               
               q=1,
               measure=list(intensity="Lamb",df=list("dnorm(z, 0, 1)")),
               measure.type="CP") 

true.parm1 <-list(a1=1.39631, a2=0.05029,
                  b0=1,b1=2,
                  Lamb=1)

# We generate a sample path.

samp1<-setSampling(Terminal=100,n=200)
set.seed(123)
sim1<-simulate(mod1,
               true.parameter=true.parm1,
               sampling=samp1)

# We estimate the parameter using qmle.
carmaopt1 <- qmle(sim1, start=true.parm1)
summary(carmaopt1)
# Internally qmle uses CarmaNoise. The result is in 
plot(carmaopt1)

# Ex.3: Carma(p=2,q=1) with scale and location parameters 
# driven by a Compound Poisson
# with jump size normally distributed.
mod2<-setCarma(p=2,                
               q=1,
               loc.par="mu",
               scale.par="sig",
               measure=list(intensity="Lamb",df=list("dnorm(z, 0, 1)")),
               measure.type="CP") 

true.parm2 <-list(a1=1.39631,
                  a2=0.05029,
                  b0=1,
                  b1=2,
                  Lamb=1,
                  mu=0.5,
                  sig=0.23)
# We simulate the sample path 
set.seed(123)
sim2<-simulate(mod2,
               true.parameter=true.parm2,
               sampling=samp1)

# We estimate the Carma and we plot the underlying noise.

carmaopt2 <- qmle(sim2, start=true.parm2)
summary(carmaopt2)

# Increments estimated by CarmaNoise
plot(carmaopt2)

## End(Not run)

Description

This function estimates the covariance between two Ito processes when they are observed at discrete times possibly nonsynchronously. It can apply to irregularly sampled one-dimensional data as a special case.

Usage

cce(x, method="HY", theta, kn, g=function(x)min(x,1-x), refreshing = TRUE,
    cwise = TRUE, delta = 0, adj = TRUE, K, c.two, J = 1, c.multi, kernel, H,
    c.RK, eta = 3/5, m = 2, ftregion = 0, vol.init = NA,
    covol.init = NA, nvar.init = NA, ncov.init = NA, mn, alpha = 0.4,
    frequency = 300, avg = TRUE, threshold, utime, psd = FALSE)

Arguments

Details

This function is a method for objects of yuima.data-class and yuima-class. It extracts the data slot when applied to a an object of yuima-class.

Typical usages are

cce(x,psd=FALSE)
cce(x,method="PHY",theta,kn,g,refreshing=TRUE,cwise=TRUE,psd=FALSE)
cce(x,method="MRC",theta,kn,g,delta=0,avg=TRUE,psd=FALSE)
cce(x,method="TSCV",K,c.two,J=1,adj=TRUE,utime,psd=FALSE)
cce(x,method="GME",c.multi,utime,psd=FALSE)
cce(x,method="RK",kernel,H,c.RK,eta=3/5,m=2,ftregion=0,utime,psd=FALSE)
cce(x,method="QMLE",vol.init=NULL,covol.init=NULL,
    nvar.init=NULL,ncov.init=NULL,psd=FALSE)
cce(x,method="SIML",mn,alpha=0.4,psd=FALSE)
cce(x,method="THY",threshold,psd=FALSE)
cce(x,method="PTHY",theta,kn,g,threshold,refreshing=TRUE,cwise=TRUE,psd=FALSE)
cce(x,method="SRC",frequency=300,avg=TRUE,utime,psd=FALSE)
cce(x,method="SBPC",frequency=300,avg=TRUE,utime,psd=FALSE)

The default method is method "HY", which is an implementation of the Hayashi-Yoshida estimator proposed in Hayashi and Yoshida (2005).

Method "PHY" is an implementation of the Pre-averaged Hayashi-Yoshida estimator proposed in Christensen et al. (2010).

Method "MRC" is an implementation of the Modulated Realized Covariance based on refresh time sampling proposed in Christensen et al. (2010).

Method "TSCV" is an implementation of the previous tick Two Scales realized CoVariance based on refresh time sampling proposed in Zhang (2011).

Method "GME" is an implementation of the Generalized Multiscale Estimator proposed in Bibinger (2011).

Method "RK" is an implementation of the multivariate Realized Kernel based on refresh time sampling proposed in Barndorff-Nielsen et al. (2011).

Method "QMLE" is an implementation of the nonparametric Quasi Maximum Likelihood Estimator proposed in Ait-Sahalia et al. (2010).

Method "SIML" is an implementation of the Separating Information Maximum Likelihood estimator proposed in Kunitomo and Sato (2013) with the basis of refresh time sampling.

Method "THY" is an implementation of the Truncated Hayashi-Yoshida estimator proposed in Mancini and Gobbi (2012).

Method "PTHY" is an implementation of the Pre-averaged Truncated Hayashi-Yoshida estimator, which is a thresholding version of the pre-averaged Hayashi-Yoshida estimator.

Method "SRC" is an implementation of the calendar time Subsampled Realized Covariance.

Method "SBPC" is an implementation of the calendar time Subsampled realized BiPower Covariation.

The rough estimation procedures for selecting the default values of the tuning parameters are based on those in Barndorff-Nielsen et al. (2009).

For the methods "PHY" or "PTHY", the default value of kn changes depending on the values of refreshing and cwise. If both refreshing and cwise are TRUE (the default), the default value of kn is given by the matrix ceiling(theta*N), where N is a matrix whose diagonal components are identical with the vector length(x)-1 and whose $(i,j)$-th component is identical with the number of the refresh times associated with $i$-th and $j$-th components of x minus 1. If refreshing is TRUE while cwise is FALSE, the default value of kn is given by ceiling(mean(theta)*sqrt(n)), where n is the number of the refresh times associated with the data minus 1. If refreshing is FALSE while cwise is TRUE, the default value of kn is given by the matrix ceiling(theta*N0), where N0 is a matrix whose diagonal components are identical with the vector length(x)-1 and whose $(i,j)$-th component is identical with (length(x)[i]-1)+(length(x)[j]-1). If both refreshing and cwise are FALSE, the default value of kn is given by ceiling(mean(theta)*sqrt(sum(length(x)-1))) (following Christensen et al. (2013)).

For the method "QMLE", the optimization of the quasi-likelihood function is implemented via arima0 using the fact that it can be seen as an MA(1) model's one: See Hansen et al. (2008) for details.

Value

A list with components:

Note

The example shows the central limit theorem for the nonsynchronous covariance estimator. Estimation of the asymptotic variance can be implemented by hyavar. The second-order correction will be provided in a future version of the package.

Author(s)

Yuta Koike with YUIMA Project Team

References

Ait-Sahalia, Y., Fan, J. and Xiu, D. (2010) High-frequency covariance estimates with noisy and asynchronous financial data, Journal of the American Statistical Association, 105, no. 492, 1504–1517.

Barndorff-Nielsen, O. E., Hansen, P. R., Lunde, A. and Shephard, N. (2008) Designing realised kernels to measure the ex-post variation of equity prices in the presence of noise, Econometrica, 76, no. 6, 1481–1536.

Barndorff-Nielsen, O. E., Hansen, P. R., Lunde, A. and Shephard, N. (2009) Realized kernels in practice: trades and quotes, Econometrics Journal, 12, C1–C32.

Barndorff-Nielsen, O. E., Hansen, P. R., Lunde, A. and Shephard, N. (2011) Multivariate realised kernels: Consistent positive semi-definite estimators of the covariation of equity prices with noise and non-synchronous trading, Journal of Econometrics, 162, 149–169.

Bibinger, M. (2011) Efficient covariance estimation for asynchronous noisy high-frequency data, Scandinavian Journal of Statistics, 38, 23–45.

Bibinger, M. (2012) An estimator for the quadratic covariation of asynchronously observed Ito processes with noise: asymptotic distribution theory, Stochastic processes and their applications, 122, 2411–2453.

Christensen, K., Kinnebrock, S. and Podolskij, M. (2010) Pre-averaging estimators of the ex-post covariance matrix in noisy diffusion models with non-synchronous data, Journal of Econometrics, 159, 116–133.

Christensen, K., Podolskij, M. and Vetter, M. (2013) On covariation estimation for multivariate continuous Ito semimartingales with noise in non-synchronous observation schemes, Journal of Multivariate Analysis 120 59–84.

Hansen, P. R., Large, J. and Lunde, A. (2008) Moving average-based estimators of integrated variance, Econometric Reviews, 27, 79–111.

Hayashi, T. and Yoshida, N. (2005) On covariance estimation of non-synchronously observed diffusion processes, Bernoulli, 11, no. 2, 359–379.

Hayashi, T. and Yoshida, N. (2008) Asymptotic normality of a covariance estimator for nonsynchronously observed diffusion processes, Annals of the Institute of Statistical Mathematics, 60, no. 2, 367–406.

Koike, Y. (2016) Estimation of integrated covariances in the simultaneous presence of nonsynchronicity, microstructure noise and jumps, Econometric Theory, 32, 533–611.

Koike, Y. (2014) An estimator for the cumulative co-volatility of asynchronously observed semimartingales with jumps, Scandinavian Journal of Statistics, 41, 460–481.

Kunitomo, N. and Sato, S. (2013) Separating information maximum likelihood estimation of realized volatility and covariance with micro-market noise, North American Journal of Economics and Finance, 26, 282–309.

Mancini, C. and Gobbi, F. (2012) Identifying the Brownian covariation from the co-jumps given discrete observations, Econometric Theory, 28, 249–273.

Varneskov, R. T. (2016) Flat-top realized kernel estimation of quadratic covariation with non-synchronous and noisy asset prices, Journal of Business & Economic Statistics, 34, no.1, 1–22.

Zhang, L. (2006) Efficient estimation of stochastic volatility using noisy observations: a multi-scale approach, Bernoulli, 12, no.6, 1019–1043.

Zhang, L. (2011) Estimating covariation: Epps effect, microstructure noise, Journal of Econometrics, 160, 33–47.

Zhang, L., Mykland, P. A. and Ait-Sahalia, Y. (2005) A tale of two time scales: Determining integrated volatility with noisy high-frequency data, Journal of the American Statistical Association, 100, no. 472, 1394–1411.

See Also

setModel, setData, hyavar, lmm, cce.factor

Examples

## Not run: 
## Set a model
diff.coef.1 <- function(t, x1 = 0, x2 = 0) sqrt(1+t)
diff.coef.2 <- function(t, x1 = 0, x2 = 0) sqrt(1+t^2)
cor.rho <- function(t, x1 = 0, x2 = 0) sqrt(1/2)
diff.coef.matrix <- matrix(c("diff.coef.1(t,x1,x2)", 
"diff.coef.2(t,x1,x2) * cor.rho(t,x1,x2)", 
"", "diff.coef.2(t,x1,x2) * sqrt(1-cor.rho(t,x1,x2)^2)"), 2, 2) 
cor.mod <- setModel(drift = c("", ""), 
diffusion = diff.coef.matrix,solve.variable = c("x1", "x2")) 

set.seed(111) 

## We use a function poisson.random.sampling to get observation by Poisson sampling.
yuima.samp <- setSampling(Terminal = 1, n = 1200) 
yuima <- setYuima(model = cor.mod, sampling = yuima.samp) 
yuima <- simulate(yuima) 
psample<- poisson.random.sampling(yuima, rate = c(0.2,0.3), n = 1000) 

## cce takes the psample and returns an estimate of the quadratic covariation. 
cce(psample)$covmat[1, 2]
##cce(psample)[1, 2]

## True value of the quadratic covariation.
cc.theta <- function(T, sigma1, sigma2, rho) { 
  tmp <- function(t) return(sigma1(t) * sigma2(t) * rho(t)) 
	integrate(tmp, 0, T) 
}

theta <- cc.theta(T = 1, diff.coef.1, diff.coef.2, cor.rho)$value 
cat(sprintf("theta =%.5f\n", theta))

names(psample@zoo.data)






# Example. A stochastic differential equation with nonlinear feedback. 

## Set a model
drift.coef.1 <- function(x1,x2) x2
drift.coef.2 <- function(x1,x2) -x1
drift.coef.vector <- c("drift.coef.1","drift.coef.2")
diff.coef.1 <- function(t,x1,x2) sqrt(abs(x1))*sqrt(1+t)
diff.coef.2 <- function(t,x1,x2) sqrt(abs(x2))
cor.rho <- function(t,x1,x2) 1/(1+x1^2)
diff.coef.matrix <- matrix(c("diff.coef.1(t,x1,x2)", 
"diff.coef.2(t,x1,x2) * cor.rho(t,x1,x2)","", 
"diff.coef.2(t,x1,x2) * sqrt(1-cor.rho(t,x1,x2)^2)"), 2, 2) 
cor.mod <- setModel(drift = drift.coef.vector,
 diffusion = diff.coef.matrix,solve.variable = c("x1", "x2"))

## Generate a path of the process
set.seed(111) 
yuima.samp <- setSampling(Terminal = 1, n = 10000) 
yuima <- setYuima(model = cor.mod, sampling = yuima.samp) 
yuima <- simulate(yuima, xinit=c(2,3)) 
plot(yuima)


## The "true" value of the quadratic covariation.
cce(yuima)

## We use the function poisson.random.sampling to generate nonsynchronous 
## observations by Poisson sampling.
psample<- poisson.random.sampling(yuima, rate = c(0.2,0.3), n = 3000) 

## cce takes the psample to return an estimated value  of the quadratic covariation. 
## The off-diagonal elements are the value of the Hayashi-Yoshida estimator. 
cce(psample)




# Example. Epps effect for the realized covariance estimator

## Set a model
drift <- c(0,0)

sigma1 <- 1
sigma2 <- 1
rho <- 0.5

diffusion <- matrix(c(sigma1,sigma2*rho,0,sigma2*sqrt(1-rho^2)),2,2)

model <- setModel(drift=drift,diffusion=diffusion,
                  state.variable=c("x1","x2"),solve.variable=c("x1","x2"))
                  
## Generate a path of the latent process
set.seed(116)

## We regard the unit interval as 6.5 hours and generate the path on it 
## with the step size equal to 2 seconds

yuima.samp <- setSampling(Terminal = 1, n = 11700) 
yuima <- setYuima(model = model, sampling = yuima.samp) 
yuima <- simulate(yuima)

## We extract nonsynchronous observations from the path generated above 
## by Poisson random sampling with the average duration equal to 10 seconds

psample <- poisson.random.sampling(yuima, rate = c(1/5,1/5), n = 11700)

## Hayashi-Yoshida estimator consistetly estimates the true correlation
cce(psample)$cormat[1,2]

## If we synchronize the observation data on some regular grid 
## by previous-tick interpolations and compute the correlation 
## by therealized covariance based on such synchronized observations, 
## we underestimate the true correlation (known as the Epps effect). 
## This is illustrated by the following examples.

## Synchronization on the grid with 5 seconds steps
suppressWarnings(s1 <- cce(subsampling(psample, sampling = setSampling(n = 4680)))$cormat[1,2])
s1

## Synchronization on the grid with 10 seconds steps
suppressWarnings(s2 <- cce(subsampling(psample, sampling = setSampling(n = 2340)))$cormat[1,2])
s2

## Synchronization on the grid with 20 seconds steps
suppressWarnings(s3 <- cce(subsampling(psample, sampling = setSampling(n = 1170)))$cormat[1,2])
s3

## Synchronization on the grid with 30 seconds steps
suppressWarnings(s4 <- cce(subsampling(psample, sampling = setSampling(n = 780)))$cormat[1,2])
s4

## Synchronization on the grid with 1 minute steps
suppressWarnings(s5 <- cce(subsampling(psample, sampling = setSampling(n = 390)))$cormat[1,2])
s5

plot(zoo(c(s1,s2,s3,s4,s5),c(5,10,20,30,60)),type="b",xlab="seconds",ylab="correlation",
main = "Epps effect for the realized covariance")



# Example. Non-synchronous and noisy observations of a correlated bivariate Brownian motion

## Generate noisy observations from the model used in the previous example
Omega <- 0.005*matrix(c(1,rho,rho,1),2,2) # covariance matrix of noise
noisy.psample <- noisy.sampling(psample,var.adj=Omega)
plot(noisy.psample)

## Hayashi-Yoshida estimator: inconsistent
cce(noisy.psample)$covmat

## Pre-averaged Hayashi-Yoshida estimator: consistent
cce(noisy.psample,method="PHY")$covmat

## Generalized multiscale estimator: consistent
cce(noisy.psample,method="GME")$covmat

## Multivariate realized kernel: consistent
cce(noisy.psample,method="RK")$covmat

## Nonparametric QMLE: consistent
cce(noisy.psample,method="QMLE")$covmat

## End(Not run)

Description

This function estimates the covariance and precision matrices of a high-dimesnional Ito process by factor modeling and regularization when it is observed at discrete times possibly nonsynchronously with noise.

Usage

cce.factor(yuima, method = "HY", factor = NULL, PCA = FALSE, 
           nfactor = "interactive", regularize = "glasso", taper, 
           group = 1:(dim(yuima) - length(factor)), lambda = "bic", 
           weight = TRUE, nlambda = 10, ratio, N, thr.type = "soft", 
           thr = NULL, tau = NULL, par.alasso = 1, par.scad = 3.7, 
           thr.delta = 0.01, frequency = 300, utime, ...)

Arguments

Details

One basic approach to estimate the covariance matrix of high-dimensional time series is to take account of the factor structure and perform regularization for the residual covariance matrix. This function implements such an estimation procedure for high-frequency data modeled as a discretely observed semimartingale. Specifically, let $Y$ be a $d$-dimensional semimartingale which describes the dynamics of the observation data. We consider the following continuous-time factor model:

$Y_t = \beta X_t + Z_t, 0\le t\le T,$

where $X$ is an $r$-dimensional semimartingale (the factor process), $Z$ is a $d$-dimensional semimartingale (the residual process), and $\beta$ is a constant $d\times r$ matrix (the factor loading matrix). We assume that $X$ and $Z$ are orthogonal in the sense that $[X,Z]_T=0$. Then, the quadratic covariation matrix of $Y$ is given by

$[Y,Y]_T=\beta[X,X]_T\beta^\top+[Z,Z]_T.$

Also, $\beta$ can be written as $\beta=[Y,X]_T[X,X]_T^{-1}$. Thus, if we have observation data both for $Y$ and $X$, we can construct estimators for $[Y,Y]_T$, $[X,X]_T$ and $\beta$ by cce. Moreover, plugging these estimators into the above equation, we can also construct an estimator for $[Z,Z]_T$. Since this estimator is often poor due to the high-dimensionality, we regularize it by some method. Then, by plugging the regularized estimator for $[Z,Z]_T$ into the above equation, we obtain the final estimator for $[Y,Y]_T$.

Even if we do not have observation data for $X$, we can (at least formally) construct a pseudo factor process by performing principal component analysis for the initial estimator of $[Y,Y]_T$. See Ait-Sahalia and Xiu (2017) and Dai et al. (2019) for details.

Currently, the following four options are available for the regularization method applied to the residual covariance matrix estimate:

1. regularize = "glasso" (the default).

   This performs the glaphical Lasso. When weight=TRUE (the default), a weighted version of the graphical Lasso is performed as in Koike (2020). Otherwise, the standard graphical Lasso is performed as in Brownlees et al. (2018).

   If lambda="aic" (resp.~lambda="bic"), the penalty parameter for the graphical Lasso is selected by minimizing the formally defined AIC (resp.~BIC). The minimization is carried out by grid search, where the grid is determined as in Section 5.1 of Koike (2020).

   The optimization problem in the graphical Lasso is solved by the GLASSOFAST algorithm of Sustik and Calderhead (2012), which is available from the package glassoFast.

2. regularize = "tapering".

   This performs tapering, i.e. taking the entry-wise product of the residual covariance matrix estimate and a tapering matrix specified by taper. See Section 3.5.1 of Pourahmadi (2011) for an overview of this method.

   If taper is missing, it is constructed according to group as follows: taper is a 0-1 matrix and the $(i,j)$-th entry is equal to 1 if and only if group[i]==group[j]. Thus, by default it makes the residual covariance matrix diagonal.

3. regularize = "thresholding".

   This performs thresholding, i.e. entries of the residual covariance matrix are shrunk toward 0 according to a thresholding rule (specified by thr.type) and a threshold level (spencified by thr).

   If thr=NULL, the $(i,j)$-th entry of thr is given by $\tau\sqrt{\hat{[Z^i,Z^i]}_T\hat{[Z^j,Z^j]}_T}$, where $\hat{[Z^i,Z^i]}_T$ (resp. $\hat{[Z^j,Z^j]}_T$) denotes the $i$-th (resp. $j$-th) diagonal entry of the non-regularized estimator for the residual covariance matrix $[Z,Z]_T$, and $\tau$ is a tuning parameter specified by tau.

   When tau=NULL, the value of $\tau$ is set to the smallest value in the grid with step size thr.delta such that the regularized estimate of the residual covariance matrix becomes positive definite.

4. regularize = "eigen.cleaning".

   This performs the eigenvalue cleaning algorithm described in Hautsch et al. (2012).

Value

A list with components:

Author(s)

Yuta Koike with YUIMA project Team

References

Ait-Sahalia, Y. and Xiu, D. (2017). Using principal component analysis to estimate a high dimensional factor model with high-frequency data, Journal of Econometrics, 201, 384–399.

Brownlees, C., Nualart, E. and Sun, Y. (2018). Realized networks, Journal of Applied Econometrics, 33, 986–1006.

Dai, C., Lu, K. and Xiu, D. (2019). Knowing factors or factor loadings, or neither? Evaluating estimators of large covariance matrices with noisy and asynchronous data, Journal of Econometrics, 208, 43–79.

Hautsch, N., Kyj, L. M. and Oomen, R. C. (2012). A blocking and regularization approach to high-dimensional realized covariance estimation, Journal of Applied Econometrics, 27, 625–645.

Koike, Y. (2020). De-biased graphical Lasso for high-frequency data, Entropy, 22, 456.

Pourahmadi, M. (2011). Covariance estimation: The GLM and regularization perspectives. Statistical Science, 26, 369–387.

Sustik, M. A. and Calderhead, B. (2012). GLASSOFAST: An efficient GLASSO implementation, UTCSTechnical Report TR-12-29, The University of Texas at Austin.

See Also

cce, lmm, glassoFast

Examples

## Not run: 
set.seed(123)

## Simulating a factor process (Heston model)
drift <- c("mu*S", "-theta*(V-v)")
diffusion <- matrix(c("sqrt(max(V,0))*S", "gamma*sqrt(max(V,0))*rho",
                       0, "gamma*sqrt(max(V,0))*sqrt(1-rho^2)"),
                    2,2)
mod <- setModel(drift = drift, diffusion = diffusion,
                state.variable = c("S", "V")) 
n <- 2340
samp <- setSampling(n = n)
heston <- setYuima(model = mod, sampling = samp)
param <- list(mu = 0.03, theta = 3, v = 0.09, 
              gamma = 0.3, rho = -0.6) 
result <- simulate(heston, xinit = c(1, 0.1), 
                   true.parameter = param)

zdata <- get.zoo.data(result) # extract the zoo data
X <- log(zdata[[1]]) # log-price process
V <- zdata[[2]] # squared volatility process


## Simulating a residual process (correlated BM)
d <- 100 # dimension
Q <- 0.1 * toeplitz(0.7^(1:d-1)) # residual covariance matrix
dZ <- matrix(rnorm(n*d),n,d) %*% chol(Q)/sqrt(n)
Z <- zoo(apply(dZ, 2, "diffinv"), samp@grid[[1]])


## Constructing observation data
b <- runif(d, 0.25, 2.25) # factor loadings
Y <- X %o% b + Z
yuima <- setData(cbind(X, Y))

# We subsample yuima to construct observation data
yuima <- subsampling(yuima, setSampling(n = 78))


## Estimating the covariance matrix (factor is known)
cmat <- tcrossprod(b) * mean(V[-1]) + Q # true covariance matrix
pmat <- solve(cmat) # true precision matrix

# (1) Regularization method is glasso (the default)
est <- cce.factor(yuima, factor = 1) 
norm(est$covmat.y - cmat, type = "2") 
norm(est$premat.y - pmat, type = "2")

# (2) Regularization method is tapering
est <- cce.factor(yuima, factor = 1, regularize = "tapering") 
norm(est$covmat.y - cmat, type = "2") 
norm(est$premat.y - pmat, type = "2")

# (3) Regularization method is thresholding
est <- cce.factor(yuima, factor = 1, regularize = "thresholding") 
norm(est$covmat.y - cmat, type = "2") 
norm(est$premat.y - pmat, type = "2")

# (4) Regularization method is eigen.cleaning
est <- cce.factor(yuima, factor = 1, regularize = "eigen.cleaning") 
norm(est$covmat.y - cmat, type = "2") 
norm(est$premat.y - pmat, type = "2")


## Estimating the covariance matrix (factor is unknown)
yuima2 <- setData(Y)

# We subsample yuima to construct observation data
yuima2 <- subsampling(yuima2, setSampling(n = 78))

# (A) Ignoring the factor structure (regularize = "glasso")
est <- cce.factor(yuima2) 
norm(est$covmat.y - cmat, type = "2") 
norm(est$premat.y - pmat, type = "2")

# (B) Estimating the factor by PCA (regularize = "glasso")
est <- cce.factor(yuima2, PCA = TRUE, nfactor = 1) # use 1 factor 
norm(est$covmat.y - cmat, type = "2") 
norm(est$premat.y - pmat, type = "2")

# One can interactively select the number of factors
# after implementing PCA (the scree plot is depicted)
# Try: est <- cce.factor(yuima2, PCA = TRUE)

## End(Not run)

Description

The yuima.PPR.qmle class is a class of the yuima package that extends the mle-class of the stats4 package.

Slots

call:

is an object of class language.

coef:

is an object of class numeric that contains estimated parameters.

fullcoef:

is an object of class numeric that contains estimated and fixed parameters.

vcov:

is an object of class matrix.

min:

is an object of class numeric.

minuslogl:

is an object of class function.

method:

is an object of class character.

model:

is an object of class yuima.PPR-class.

Methods

Methods mle

All methods for mle-class are available.

Author(s)

The YUIMA Project Team

Description

The cogarch.est class is a class of the yuima package that contains estimated parameters obtained by the function gmm or qmle.

Slots

yuima:

is an object of of yuima-class.

objFun:

is an object of class character that indicates the objective function used in the minimization problem. See the documentation of the function gmm or qmle for more details.

call:

is an object of class language.

coef:

is an object of class numeric that contains estimated parameters.

fullcoef:

is an object of class numeric that contains estimated and fixed parameters.

vcov:

is an object of class matrix.

min:

is an object of class numeric.

minuslogl:

is an object of class function.

method:

is an object of class character.

Methods

Methods mle

All methods for mle-class are available.

Author(s)

The YUIMA Project Team

Description

The cogarch.est.incr class is a class of the yuima package that extends the cogarch.est-class and is filled by the function gmm or qmle.

Slots

Incr.Lev:

is an object of class zoo that contains the estimated increments of the noise obtained using cogarchNoise.

yuima:

is an object of of yuima-class.

logL.Incr:

is an object of class numeric that contains the value of the log-likelihood for estimated Levy increments.

objFun:

is an object of class character that indicates the objective function used in the minimization problem. See the documentation of the function gmm or qmle for more details.

call:

is an object of class language.

coef:

is an object of class numeric that contains estimated parameters.

fullcoef:

is an object of class numeric that contains estimated and fixed parameters.

vcov:

is an object of class matrix.

min:

is an object of class numeric.

minuslogl:

is an object of class function.

method:

is an object of class character.

Methods

simulate

simulation method. For more information see simulate.

plot

Plot method for estimated increment of the noise.

Methods mle

All methods for mle-class are available.

Author(s)

The YUIMA Project Team

Description

The cogarch.info-class is a class of the yuima package

Slots

p:

Number of autoregressive coefficients in the variance process.

q:

Number of moving average coefficients in the variance process.

ar.par:

Label of autoregressive coefficients.

ma.par:

Label of moving average coefficients.

loc.par:

Label of location coefficient in the variance process.

Cogarch.var:

Label of the observed process.

V.var:

Label of the variance process.

Latent.var:

Label of the latent process in the state representation of the variance.

XinExpr:

Logical variable. If XinExpr=FALSE, the starting condition of Latent.var is zero otherwise each component of Latent.var has a parameter as a starting point.

measure:

Levy measure for jump and quadratic part.

measure.type:

Type specification for Levy measure.

Note

The cogarch.info-class object cannot be directly specified by the user but it is built when the yuima.cogarch-class object is constructed via setCogarch.

Author(s)

The YUIMA Project Team

Description

Retrieve the increment of the underlying Levy for the COGARCH(p,q) process

Usage

cogarchNoise(yuima, data=NULL, param, mu=1)

Arguments

Value

Note

The function cogarchNoise assumes the underlying Levy process is centered in zero.

The function gmm uses the function cogarchNoise for estimation of underlying Levy in the COGARCH(p,q) model.

Author(s)

The YUIMA Project Team

References

Chadraa. (2009) Statistical Modelling with COGARCH(P,Q) Processes, PhD Thesis.

Examples

# Insert here some examples

Description

Volatility structural change point estimator

Usage

CPoint(yuima, param1, param2, print=FALSE, symmetrized=FALSE, plot=FALSE)
qmleL(yuima, t,  ...)
qmleR(yuima, t,  ...)

Arguments

Details

CPoint estimates the change point using quasi-maximum likelihood approach.

Function qmleL estimates the parameters in the diffusion matrix using observations up to time t.

Function qmleR estimates the parameters in the diffusion matrix using observations from time t to the end.

Arguments in both qmleL and qmleR follow the same rules as in qmle.

Value

Author(s)

The YUIMA Project Team

Examples

## Not run: 
diff.matrix <- matrix(c("theta1.1*x1","0*x2","0*x1","theta1.2*x2"), 2, 2)

drift.c <- c("1-x1", "3-x2")
drift.matrix <- matrix(drift.c, 2, 1)

ymodel <- setModel(drift=drift.matrix, diffusion=diff.matrix, time.variable="t",
state.variable=c("x1", "x2"), solve.variable=c("x1", "x2"))
n <- 1000

set.seed(123)

t1 <- list(theta1.1=.1, theta1.2=0.2)
t2 <- list(theta1.1=.6, theta1.2=.6)

tau <- 0.4
ysamp1 <- setSampling(n=tau*n, Initial=0, delta=0.01)
yuima1 <- setYuima(model=ymodel, sampling=ysamp1)
yuima1 <- simulate(yuima1, xinit=c(1, 1), true.parameter=t1)

x1 <- yuima1@data@zoo.data[[1]]
x1 <- as.numeric(x1[length(x1)])
x2 <- yuima1@data@zoo.data[[2]]
x2 <- as.numeric(x2[length(x2)])

ysamp2 <- setSampling(Initial=n*tau*0.01, n=n*(1-tau), delta=0.01)
yuima2 <- setYuima(model=ymodel, sampling=ysamp2)

yuima2 <- simulate(yuima2, xinit=c(x1, x2), true.parameter=t2)


yuima <- yuima1
yuima@data@zoo.data[[1]] <- c(yuima1@data@zoo.data[[1]], yuima2@data@zoo.data[[1]][-1])
yuima@data@zoo.data[[2]] <- c(yuima1@data@zoo.data[[2]], yuima2@data@zoo.data[[2]][-1])

plot(yuima)

# estimation of change point for given parameter values
t.est <- CPoint(yuima,param1=t1,param2=t2, plot=TRUE)


low <- list(theta1.1=0, theta1.2=0)

# first state estimate of parameters using small 
# portion of data in the tails
tmp1 <- qmleL(yuima,start=list(theta1.1=0.3,theta1.2=0.5),t=1.5,
        lower=low, method="L-BFGS-B")
tmp1
tmp2 <- qmleR(yuima,start=list(theta1.1=0.3,theta1.2=0.5), t=8.5,
        lower=low, method="L-BFGS-B")
tmp2


# first stage changepoint estimator
t.est2 <- CPoint(yuima,param1=coef(tmp1),param2=coef(tmp2))
t.est2$tau


# second stage estimation of parameters given first stage
# change point estimator
tmp11 <- qmleL(yuima,start=as.list(coef(tmp1)), t=t.est2$tau-0.1,
 lower=low, method="L-BFGS-B")
tmp11

tmp21 <- qmleR(yuima,start=as.list(coef(tmp2)), t=t.est2$tau+0.1,
 lower=low, method="L-BFGS-B")
tmp21

# second stage estimator of the change point
CPoint(yuima,param1=coef(tmp11),param2=coef(tmp21))


## One dimensional example: non linear case
diff.matrix <- matrix("(1+x1^2)^theta1", 1, 1)
drift.c <- c("x1")

ymodel <- setModel(drift=drift.c, diffusion=diff.matrix, time.variable="t",
state.variable=c("x1"), solve.variable=c("x1"))
n <- 500

set.seed(123)

y0 <- 5 # initial value
theta00 <- 1/5
gamma <- 1/4


theta01 <- theta00+n^(-gamma)


t1 <- list(theta1= theta00)
t2 <- list(theta1= theta01)

tau <- 0.4
ysamp1 <- setSampling(n=tau*n, Initial=0, delta=1/n)
yuima1 <- setYuima(model=ymodel, sampling=ysamp1)
yuima1 <- simulate(yuima1, xinit=c(5), true.parameter=t1)
x1 <- yuima1@data@zoo.data[[1]]
x1 <- as.numeric(x1[length(x1)])

ysamp2 <- setSampling(Initial=tau, n=n*(1-tau), delta=1/n)
yuima2 <- setYuima(model=ymodel, sampling=ysamp2)

yuima2 <- simulate(yuima2, xinit=c(x1), true.parameter=t2)


yuima <- yuima1
yuima@data@zoo.data[[1]] <- c(yuima1@data@zoo.data[[1]], yuima2@data@zoo.data[[1]][-1])


plot(yuima)


t.est <- CPoint(yuima,param1=t1,param2=t2)
t.est$tau

low <- list(theta1=0)
upp <- list(theta1=1)

# first state estimate of parameters using small 
# portion of data in the tails
tmp1 <- qmleL(yuima,start=list(theta1=0.5),t=.15,lower=low, upper=upp,method="L-BFGS-B")
tmp1
tmp2 <- qmleR(yuima,start=list(theta1=0.5), t=.85,lower=low, upper=upp,method="L-BFGS-B")
tmp2



# first stage changepoint estimator
t.est2 <- CPoint(yuima,param1=coef(tmp1),param2=coef(tmp2))
t.est2$tau


# second stage estimation of parameters given first stage
# change point estimator
tmp11 <- qmleL(yuima,start=as.list(coef(tmp1)), t=t.est2$tau-0.1,
   lower=low, upper=upp,method="L-BFGS-B")
tmp11

tmp21 <- qmleR(yuima,start=as.list(coef(tmp2)), t=t.est2$tau+0.1,
  lower=low, upper=upp,method="L-BFGS-B")
tmp21


# second stage estimator of the change point
CPoint(yuima,param1=coef(tmp11),param2=coef(tmp21),plot=TRUE)

## End(Not run)

Description

The function converts an object of class zoo to an object of class yuima.PPR.

Usage

DataPPR(CountVar, yuimaPPR, samp)

Arguments

Value

The function returns an object of class yuima.PPR where the slot model contains the Point process described in yuimaPPR@model, the slot data contains the counting variables and the covariates observed on the grid in samp.

Examples

## Not run: 
# In this example we generate a dataset contains the Counting Variable N
# and the Covariate X.
# The covariate X is an OU driven by a Gamma process.

# Values of parameters.
mu <- 2
alpha <- 4
beta <-5

# Law definition
my.rKern <- function(n,t){
  res0 <- t(t(rgamma(n, 0.1*t)))
  res1 <- t(t(rep(1,n)))
  res <- cbind(res0,res1)
  return(res)
}

Law.PPRKern <- setLaw(rng = my.rKern)

# Point Process definition
modKern <- setModel(drift = c("0.4*(0.1-X)","0"),
                    diffusion = c("0","0"),
                    jump.coeff = matrix(c("1","0","0","1"),2,2),
                    measure = list(df = Law.PPRKern),
                    measure.type = c("code","code"),
                    solve.variable = c("X","N"),
                    xinit=c("0.25","0"))

gFun <- "exp(mu*log(1+X))"
#
Kernel <- "alpha*exp(-beta*(t-s))"

prvKern <- setPPR(yuima = modKern,
                  counting.var="N", gFun=gFun,
                  Kernel = as.matrix(Kernel),
                  lambda.var = "lambda", var.dx = "N",
                  lower.var="0", upper.var = "t")

# Simulation

Term<-200
seed<-1
n<-20000

true.parKern <- list(mu=mu, alpha=alpha, beta=beta)


set.seed(seed)
# set.seed(1)

time.simKern <-system.time(
  simprvKern <- simulate(object = prvKern, true.parameter = true.parKern,
                         sampling = setSampling(Terminal =Term, n=n))
)


plot(simprvKern,main ="Counting Process with covariates" ,cex.main=0.9)

# Using the function get.counting.data we extract from an object of class
# yuima.PPR the counting process N and the covariate X at the arrival times.

CountVar <- get.counting.data(simprvKern)

plot(CountVar)

# We convert the zoo object in the yuima.PPR object.

sim2 <- DataPPR(CountVar, yuimaPPR=simprvKern, samp=simprvKern@sampling)


## End(Not run)

Description

This function verifies if the condition of stationarity is satisfied.

Usage

Diagnostic.Carma(carma)

Arguments

Value

Logical variable. If TRUE, Carma is stationary.

Author(s)

YUIMA TEAM

Examples

mod1 <- setCarma(p = 2, q = 1, scale.par = "sig",
          Carma.var = "y")

param1 <- list(a1 = 1.39631, a2 = 0.05029, b0 = 1,
            b1 = 1, sig = 1)
samp1 <- setSampling(Terminal = 100, n = 200)

set.seed(123)

sim1 <- simulate(mod1, true.parameter = param1,
          sampling = samp1)

est1 <- qmle(sim1, start = param1)

Diagnostic.Carma(est1)

Description

The function check the statistical properties of the COGARCH(p,q) model. We verify if the process has a strict positive stationary variance model.

Usage

Diagnostic.Cogarch(yuima.cogarch, param = list(), matrixS = NULL, mu = 1, display = TRUE)

Arguments

Value

The functon returns a List with entries:

Author(s)

YUIMA Project Team

Examples

## Not run: 
# Definition of the COGARCH(1,1) process driven by a Variance Gamma nois:
param.VG <- list(a1 = 0.038,  b1 =  0.053,
                  a0 = 0.04/0.053,lambda = 1, alpha = sqrt(2), beta = 0, mu = 0, 
                  x01 = 50.33)

cog.VG <- setCogarch(p = 1, q = 1, work = FALSE,
                      measure=list(df="rvgamma(z, lambda, alpha, beta, mu)"),
                      measure.type = "code", 
                      Cogarch.var = "y",
                      V.var = "v", Latent.var="x",
                      XinExpr=TRUE)

# Verify the stationarity and the positivity of th variance process

test <- Diagnostic.Cogarch(cog.VG,param=param.VG)
show(test)

# Simulate a sample path

set.seed(210)

Term=800
num=24000

samp.VG <- setSampling(Terminal=Term, n=num)

sim.VG <- simulate(cog.VG,
                    true.parameter=param.VG,
                    sampling=samp.VG,
                    method="euler")
plot(sim.VG)

# Estimate the model

res.VG <- gmm(sim.VG, start = param.VG, Est.Incr = "IncrPar")

summary(res.VG)

#  Check if the estimated COGARCH(1,1) has a positive and stationary variance

test1<-Diagnostic.Cogarch(res.VG)
show(test1)

# Simulate a COGARCH sample path using the estimated COGARCH(1,1) 
# and the recovered increments of underlying Variance Gamma Noise

esttraj<-simulate(res.VG)
plot(esttraj)



## End(Not run)

Description

The function estimates a t-Levy Regression Model

Usage

estimation_LRM(start, model, data, upper, lower, PT = 500, n_obs1 = NULL)

Arguments

Details

A two-step estimation procedure. Regressor coefficients and scale parameters are obtained by maximizing the quasi-likelihood function based on the Cauchy density. The degree of freedom is estimated used the unitary increment of the t-noise.

Value

Estimated parameters

Author(s)

The YUIMA Project Team

Contacts: Lorenzo Mercuri [email protected]

Description

The function provides two estimation procedures: Maximum Likelihood Estimation and Matching Empirical Correlation

Usage

EstimCarmaHawkes(yuima, start, est.method = "qmle", method = "BFGS", 
lower = NULL, upper = NULL, lags = NULL, display = FALSE)

Arguments

Value

The output contains the estimated parameters.

Author(s)

The YUIMA Project Team

Contacts: Lorenzo Mercuri [email protected]

References

Mercuri, L., Perchiazzo, A., & Rroji, E. (2022). A Hawkes model with CARMA (p, q) intensity. doi:10.48550/arXiv.2208.02659.

Examples

## Not run: 
## MLE For A CARMA(2,1)-Hawkes ##

# Inputs:
a <- c(3,2)
b <- c(1,0.3)
mu<-0.30

true.par<-c(mu,a,b)

# step 1) Model Definition => Constructor 'setCarmaHawkes'
p <- 2
q <- 1
mod1 <- setCarmaHawkes(p = p,q = q)

# step 2) Grid Construction => Constructor 'setSampling'
FinalTime <- 5000
t0 <- 0
samp <- setSampling(t0, FinalTime, n = FinalTime)

# step 3) Simulation => method 'simulate'
# We use method 'simulate' to generate our dataset. 
# For the estimation from real data, 
# we use the constructors 'setData' and 
#'setYuima' (input 'model' is an object of 
#           'yuima.CarmaHawkes-class'). 

names(true.par) <- c(mod1@info@base.Int, mod1@info@ar.par, mod1@info@ma.par) 

set.seed(1)
system.time(
sim1 <- simulate(object = mod1, true.parameter = true.par, 
    sampling = samp)
)

plot(sim1)

# step 4) Estimation using the likelihood function.
system.time(
  res <- EstimCarmaHawkes(yuima = sim1, 
    start = true.par)
)


## End(Not run)

Description

This is a function to simulate the preliminary estimator and the corresponding one step estimators based on the Newton-Raphson and the scoring method of the Cox-Ingersoll-Ross process given via the SDE

$\mathrm{d} X_t = (\alpha-\beta X_t)\mathrm{d} t + \sqrt{\gamma X_t}\mathrm{d} W_t$

with parameters $\beta>0,$ $2\alpha>5\gamma>0$ and a Brownian motion $(W_t)_{t\geq 0}$. This function uses the Gaussian quasi-likelihood, hence requires that data is sampled at high-frequency.

Usage

fitCIR(data)

Arguments

Details

The estimators calculated by this function can be found in the reference below.

Value

A list with three entries each contain a vector in the following order: The result of the preliminary estimator, Newton-Raphson method and the method of scoring.

If the sampling points are not equidistant the function will return 'Please use equidistant sampling points'.

Author(s)

Nicole Hufnagel

Contacts: [email protected]

References

Y. Cheng, N. Hufnagel, H. Masuda. Estimation of ergodic square-root diffusion under high-frequency sampling. Econometrics and Statistics, Article Number: 346 (2022).

Examples

#You can make use of the function simCIR to generate the data 
data <- simCIR(alpha=3,beta=1,gamma=1, n=5000, h=0.05, equi.dist=TRUE)
results <- fitCIR(data)