[Yuima-commits] r109 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jul 20 13:11:25 CEST 2010
Author: iacus
Date: 2010-07-20 13:11:25 +0200 (Tue, 20 Jul 2010)
New Revision: 109
Added:
pkg/yuima/R/CPoint.R
pkg/yuima/man/CPoint.Rd
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/NEWS
pkg/yuima/R/qmle.R
Log:
added CPoint
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2010-07-19 11:18:15 UTC (rev 108)
+++ pkg/yuima/DESCRIPTION 2010-07-20 11:11:25 UTC (rev 109)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package (unstable version)
-Version: 0.1.06
-Date: 2010-07-19
+Version: 0.1.07
+Date: 2010-07-20
Depends: methods, zoo, stats4, utils
Suggests: cubature, mvtnorm
Author: YUIMA Project Team.
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2010-07-19 11:18:15 UTC (rev 108)
+++ pkg/yuima/NAMESPACE 2010-07-20 11:11:25 UTC (rev 109)
@@ -91,6 +91,9 @@
export(qmle)
export(quasilogl)
export(lasso)
+export(CPoint)
+export(qmleR)
+export(qmleL)
S3method(toLatex, yuima)
Modified: pkg/yuima/NEWS
===================================================================
--- pkg/yuima/NEWS 2010-07-19 11:18:15 UTC (rev 108)
+++ pkg/yuima/NEWS 2010-07-20 11:11:25 UTC (rev 109)
@@ -1,3 +1,7 @@
+Changes in Version 0.1.07
+
+ o added CPoint() for change point estimation
+
Changes in Version 0.1.05
o added lasso() for adaptive lasso estimation
Added: pkg/yuima/R/CPoint.R
===================================================================
--- pkg/yuima/R/CPoint.R (rev 0)
+++ pkg/yuima/R/CPoint.R 2010-07-20 11:11:25 UTC (rev 109)
@@ -0,0 +1,398 @@
+qmleL <- function(yuima, start, method="BFGS", t, print=FALSE, lower, upper, ...){
+
+ call <- match.call()
+
+ if( missing(yuima))
+ yuima.stop("yuima object is missing.")
+
+## param handling
+
+## FIXME: maybe we should choose initial values at random within lower/upper
+## at present, qmle stops
+ if( missing(start) )
+ yuima.stop("Starting values for the parameters are missing.")
+
+ diff.par <- yuima at model@parameter at diffusion
+
+ if(!is.list(start))
+ yuima.stop("Argument 'start' must be of list type.")
+
+ fullcoef <- diff.par
+ npar <- length(fullcoef)
+
+ nm <- names(start)
+ oo <- match(nm, fullcoef)
+ if(any(is.na(oo)))
+ yuima.stop("some named arguments in 'start' are not arguments to the supplied yuima model")
+ start <- start[order(oo)]
+ nm <- names(start)
+
+ idx.diff <- match(diff.par, nm)
+
+ tmplower <- as.list( rep( -Inf, length(nm)))
+ names(tmplower) <- nm
+ if(!missing(lower)){
+ idx <- match(names(lower), names(tmplower))
+ if(any(is.na(idx)))
+ yuima.stop("names in 'lower' do not match names fo parameters")
+ tmplower[ idx ] <- lower
+ }
+ lower <- tmplower
+
+ tmpupper <- as.list( rep( Inf, length(nm)))
+ names(tmpupper) <- nm
+ if(!missing(upper)){
+ idx <- match(names(upper), names(tmpupper))
+ if(any(is.na(idx)))
+ yuima.stop("names in 'lower' do not match names fo parameters")
+ tmpupper[ idx ] <- upper
+ }
+ upper <- tmpupper
+
+ d.size <- yuima at model@equation.number
+ n <- length(yuima)[1]
+
+ times <- time(yuima at data@zoo.data[[1]])
+ minT <- as.numeric(times[1])
+ maxT <- as.numeric(times[length(times)])
+
+ if(missing(t))
+ t <- mean(c(minT,maxT))
+ k <- max(which(times <=t),na.rm=TRUE)[1]
+
+ if(k<2)
+ k <- 2
+ if(k>n-2)
+ k <- n-2
+
+ env <- new.env()
+# assign("X", as.matrix(onezoo(yuima)), env=env)
+# assign("deltaX", matrix(0, n-1, d.size), env=env)
+# for(t in 1:(n-1))
+# env$deltaX[t,] <- env$X[t+1,] - env$X[t,]
+
+ assign("X", as.matrix(onezoo(yuima))[1:k,], env=env)
+ assign("deltaX", matrix(0, k-1, d.size), env=env)
+ for(t in 1:(k-1))
+ env$deltaX[t,] <- env$X[t+1,] - env$X[t,]
+
+
+ assign("h", deltat(yuima at data@zoo.data[[1]]), env=env)
+
+ f <- function(p) {
+ mycoef <- as.list(p)
+ names(mycoef) <- nm
+ sum(pminusquasilogl(yuima=yuima, param=mycoef, print=print, env))
+ }
+
+
+ oout <- NULL
+
+
+ if(length(start)>1){ # multidimensional optim
+ oout <- optim(start, f, method = method, hessian = FALSE, lower=lower, upper=upper)
+
+ } else { ### one dimensional optim
+ opt <- optimize(f, ...) ## an interval should be provided
+ oout <- list(par = opt$minimum, value = opt$objective)
+ } ### endif( length(start)>1 )
+
+
+ oout <- oout[c("par","value")]
+ return(oout)
+}
+
+
+qmleR <- function(yuima, start, method="BFGS", t, print=FALSE, lower, upper, ...){
+
+ call <- match.call()
+
+ if( missing(yuima))
+ yuima.stop("yuima object is missing.")
+
+## param handling
+
+## FIXME: maybe we should choose initial values at random within lower/upper
+## at present, qmle stops
+ if( missing(start) )
+ yuima.stop("Starting values for the parameters are missing.")
+
+ diff.par <- yuima at model@parameter at diffusion
+
+ if(!is.list(start))
+ yuima.stop("Argument 'start' must be of list type.")
+
+ fullcoef <- diff.par
+ npar <- length(fullcoef)
+
+ nm <- names(start)
+ oo <- match(nm, fullcoef)
+ if(any(is.na(oo)))
+ yuima.stop("some named arguments in 'start' are not arguments to the supplied yuima model")
+ start <- start[order(oo)]
+ nm <- names(start)
+
+ idx.diff <- match(diff.par, nm)
+
+ tmplower <- as.list( rep( -Inf, length(nm)))
+ names(tmplower) <- nm
+ if(!missing(lower)){
+ idx <- match(names(lower), names(tmplower))
+ if(any(is.na(idx)))
+ yuima.stop("names in 'lower' do not match names fo parameters")
+ tmplower[ idx ] <- lower
+ }
+ lower <- tmplower
+
+ tmpupper <- as.list( rep( Inf, length(nm)))
+ names(tmpupper) <- nm
+ if(!missing(upper)){
+ idx <- match(names(upper), names(tmpupper))
+ if(any(is.na(idx)))
+ yuima.stop("names in 'lower' do not match names fo parameters")
+ tmpupper[ idx ] <- upper
+ }
+ upper <- tmpupper
+
+ d.size <- yuima at model@equation.number
+ n <- length(yuima)[1]
+
+ times <- time(yuima at data@zoo.data[[1]])
+ minT <- as.numeric(times[1])
+ maxT <- as.numeric(times[length(times)])
+
+ if(missing(t))
+ t <- mean(c(minT,maxT))
+ k <- max(which(times <=t),na.rm=TRUE)[1]
+
+ if(k<2)
+ k <- 2
+ if(k>n-2)
+ k <- n-2
+
+
+ env <- new.env()
+# assign("X", as.matrix(onezoo(yuima)), env=env)
+# assign("deltaX", matrix(0, n-1, d.size), env=env)
+# for(t in 1:(n-1))
+# env$deltaX[t,] <- env$X[t+1,] - env$X[t,]
+
+ assign("X", as.matrix(onezoo(yuima))[-(1:k),], env=env)
+ assign("deltaX", matrix(0, dim(env$X)[1]-1, d.size), env=env)
+ for(t in 1:(dim(env$X)[1]-1))
+ env$deltaX[t,] <- env$X[t+1,] - env$X[t,]
+
+
+ assign("h", deltat(yuima at data@zoo.data[[1]]), env=env)
+
+ f <- function(p) {
+ mycoef <- as.list(p)
+ names(mycoef) <- nm
+ sum(pminusquasilogl(yuima=yuima, param=mycoef, print=print, env))
+ }
+
+
+ oout <- NULL
+
+
+ if(length(start)>1){ # multidimensional optim
+ oout <- optim(start, f, method = method, hessian = FALSE, lower=lower, upper=upper)
+
+ } else { ### one dimensional optim
+ opt <- optimize(f, ...) ## an interval should be provided
+ oout <- list(par = opt$minimum, value = opt$objective)
+ } ### endif( length(start)>1 )
+
+
+ oout <- oout[c("par","value")]
+ return(oout)
+}
+
+CPoint <- function(yuima, param1, param2, print=FALSE, plot=FALSE){
+ d.size <- yuima at model@equation.number
+ n <- length(yuima)[1]
+
+ d.size <- yuima at model@equation.number
+
+ env <- new.env()
+ assign("X", as.matrix(onezoo(yuima)), env=env)
+ assign("deltaX", matrix(0, n-1, d.size), env=env)
+ for(t in 1:(n-1))
+ env$deltaX[t,] <- env$X[t+1,] - env$X[t,]
+
+ assign("h", deltat(yuima at data@zoo.data[[1]]), env=env)
+
+ QL1 <- pminusquasilogl(yuima=yuima, param=param1, print=print, env)
+ QL2 <- pminusquasilogl(yuima=yuima, param=param2, print=print, env)
+
+ D <- sapply(2:(n-1), function(x) sum(QL1[1:x]) + sum(QL2[-(1:x)]))
+ D <- c(D[1], D, D[length(D)])
+ D <- ts(D, start=0, deltat=deltat(yuima at data@zoo.data[[1]]))
+ if(plot)
+ plot(D,type="l", main="change point statistics")
+ tau.hat <- index(yuima at data@zoo.data[[1]])[which.min(D)]
+
+ return(list(tau=tau.hat, param1=param1, param2=param2))
+}
+
+
+
+
+# partial quasi-likelihood
+# returns a vector of conditionational minus-log-likelihood terms
+# the whole negative log likelihood is the sum
+
+pminusquasilogl <- function(yuima, param, print=FALSE, env){
+
+ diff.par <- yuima at model@parameter at diffusion
+ fullcoef <- diff.par
+ npar <- length(fullcoef)
+
+ nm <- names(param)
+ oo <- match(nm, fullcoef)
+ if(any(is.na(oo)))
+ yuima.stop("some named arguments in 'param' are not arguments to the supplied yuima model")
+
+
+ param <- param[order(oo)]
+ nm <- names(param)
+
+ idx.diff <- match(diff.par, nm)
+
+ h <- env$h
+
+ theta1 <- unlist(param[idx.diff])
+ n.theta1 <- length(theta1)
+ n.theta <- n.theta1
+
+ d.size <- yuima at model@equation.number
+#n <- length(yuima)[1]
+ n <- dim(env$X)[1]
+
+ vec <- env$deltaX
+
+ K <- -0.5*d.size * log( (2*pi*h) )
+
+
+
+ QL <- 0
+ pn <- numeric(n-1)
+ diff <- diffusion.term(yuima, param, env)
+ dimB <- dim(diff[, , 1])
+
+ if(is.null(dimB)){ # one dimensional X
+ for(t in 1:(n-1)){
+ yB <- diff[, , t]^2
+ logdet <- log(yB)
+ pn[t] <- K - 0.5*logdet-0.5*vec[t, ]^2/(h*yB)
+ QL <- QL+pn[t]
+
+ }
+ } else { # multidimensional X
+ for(t in 1:(n-1)){
+ yB <- diff[, , t] %*% t(diff[, , t])
+ logdet <- log(det(yB))
+ if(is.infinite(logdet) ){ # should we return 1e10?
+ pn[t] <- log(1)
+ yuima.warn("singular diffusion matrix")
+ return(1e10)
+ }else{
+ pn[t] <- K - 0.5*logdet +
+ ((-1/(2*h))*t(vec[t, ])%*%solve(yB)%*%vec[t, ])
+ QL <- QL+pn[t]
+ }
+ }
+ }
+
+ if(!is.finite(QL)){
+ yuima.warn("quasi likelihood is too small to calculate.")
+ QL <- 1e10
+ }
+
+ if(print==TRUE){
+ yuima.warn(sprintf("NEG-QL: %f, %s", -QL, paste(names(param),param,sep="=",collapse=", ")))
+ }
+
+
+ return(-pn)
+
+}
+
+
+pminusquasiloglL <- function(yuima, param, print=FALSE, env){
+
+ diff.par <- yuima at model@parameter at diffusion
+ fullcoef <- diff.par
+ npar <- length(fullcoef)
+
+ nm <- names(param)
+ oo <- match(nm, fullcoef)
+ if(any(is.na(oo)))
+ yuima.stop("some named arguments in 'param' are not arguments to the supplied yuima model")
+
+
+ param <- param[order(oo)]
+ nm <- names(param)
+
+ idx.diff <- match(diff.par, nm)
+
+ h <- env$h
+
+ theta1 <- unlist(param[idx.diff])
+ n.theta1 <- length(theta1)
+ n.theta <- n.theta1
+
+ d.size <- yuima at model@equation.number
+ n <- length(yuima)[1]
+
+ vec <- env$deltaX
+
+ K <- -0.5*d.size * log( (2*pi*h) )
+
+
+
+ QL <- 0
+ pn <- numeric(n-1)
+ diff <- diffusion.term(yuima, param, env)
+ dimB <- dim(diff[, , 1])
+
+ if(is.null(dimB)){ # one dimensional X
+ for(t in 1:(n-1)){
+ yB <- diff[, , t]^2
+ logdet <- log(yB)
+ pn[t] <- K - 0.5*logdet-0.5*vec[t, ]^2/(h*yB)
+ QL <- QL+pn[t]
+
+ }
+ } else { # multidimensional X
+ for(t in 1:(n-1)){
+ yB <- diff[, , t] %*% t(diff[, , t])
+ logdet <- log(det(yB))
+ if(is.infinite(logdet) ){ # should we return 1e10?
+ pn[t] <- log(1)
+ yuima.warn("singular diffusion matrix")
+ return(1e10)
+ }else{
+ pn[t] <- K - 0.5*logdet +
+ ((-1/(2*h))*t(vec[t, ])%*%solve(yB)%*%vec[t, ])
+ QL <- QL+pn[t]
+ }
+ }
+ }
+
+ if(!is.finite(QL)){
+ yuima.warn("quasi likelihood is too small to calculate.")
+ QL <- 1e10
+ }
+
+ if(print==TRUE){
+ yuima.warn(sprintf("NEG-QL: %f, %s", -QL, paste(names(param),param,sep="=",collapse=", ")))
+ }
+
+ return(-pn)
+
+}
+
+
+
+
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2010-07-19 11:18:15 UTC (rev 108)
+++ pkg/yuima/R/qmle.R 2010-07-20 11:11:25 UTC (rev 109)
@@ -15,7 +15,9 @@
d.size <- yuima at model@equation.number
modelstate <- yuima at model@state.variable
DRIFT <- yuima at model@drift
- n <- length(yuima)[1]
+# n <- length(yuima)[1]
+ n <- dim(env$X)[1]
+
drift <- matrix(0,n,d.size)
@@ -39,7 +41,9 @@
d.size <- yuima at model@equation.number
modelstate <- yuima at model@state.variable
DIFFUSION <- yuima at model@diffusion
- n <- length(yuima)[1]
+# n <- length(yuima)[1]
+ n <- dim(env$X)[1]
+
diff <- array(0, dim=c(d.size, r.size, n))
for(i in 1:length(theta)){
assign(names(theta)[i],theta[[i]])
@@ -109,10 +113,16 @@
if(!is.list(start))
yuima.stop("Argument 'start' must be of list type.")
- fullcoef <- c(diff.par, drift.par)
+ fullcoef <- NULL
+
+ if(length(diff.par)>0)
+ fullcoef <- diff.par
+
+ if(length(drift.par)>0)
+ fullcoef <- c(fullcoef, drift.par)
+
npar <- length(fullcoef)
-
fixed.par <- names(fixed)
if (any(!(fixed.par %in% fullcoef)))
@@ -167,7 +177,12 @@
f <- function(p) {
mycoef <- as.list(p)
- names(mycoef) <- nm[-idx.fixed]
+# print(nm[-idx.fixed])
+# print(nm)
+ if(length(idx.fixed)>0)
+ names(mycoef) <- nm[-idx.fixed]
+ else
+ names(mycoef) <- nm
mycoef[fixed.par] <- fixed
minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
}
@@ -198,12 +213,19 @@
oout <- list(par = opt1$minimum, value = opt1$objective)
} ### endif( length(start)>1 )
} else { ### first diffusion, then drift
+ theta1 <- NULL
+
+ old.fixed <- fixed
+ old.start <- start
+
+ if(length(idx.diff)>0){
## DIFFUSION ESTIMATIOn first
old.fixed <- fixed
old.start <- start
new.start <- start[idx.diff] # considering only initial guess for diffusion
new.fixed <- fixed
- new.fixed[nm[idx.drift]] <- start[idx.drift]
+ if(length(idx.drift)>0)
+ new.fixed[nm[idx.drift]] <- start[idx.drift]
fixed <- new.fixed
fixed.par <- names(fixed)
idx.fixed <- match(fixed.par, nm)
@@ -214,6 +236,7 @@
mydots$fn <- as.name("f")
mydots$start <- NULL
mydots$par <- unlist(new.start)
+# print(mydots)
mydots$hessian <- FALSE
mydots$upper <- unlist( upper[ nm[idx.diff] ])
mydots$lower <- unlist( lower[ nm[idx.diff] ])
@@ -234,7 +257,11 @@
oout <- list(par = theta1, value = opt1$objective)
}
theta1 <- oout$par
+ } ## endif(length(idx.diff)>0)
+ theta2 <- NULL
+
+ if(length(idx.drift)>0){
## DRIFT estimation with first state diffusion estimates
fixed <- old.fixed
start <- old.start
@@ -270,7 +297,8 @@
oout1 <- list(par = theta2, value = as.numeric(opt1$objective))
}
theta2 <- oout1$par
- oout1$par <- c(theta1, theta2)
+ } ## endif(length(idx.drift)>0)
+ oout1 <- list(par= c(theta1, theta2))
names(oout1$par) <- c(diff.par,drift.par)
oout <- oout1
@@ -300,6 +328,8 @@
names(par) <- c(diff.par, drift.par)
nm <- c(diff.par, drift.par)
+# print(par)
+# print(coef)
conDrift <- list(trace = 5, fnscale = 1,
parscale = rep.int(5, length(drift.par)),
ndeps = rep.int(0.001, length(drift.par)), maxit = 100L,
@@ -335,12 +365,12 @@
# if (npar == 1 && method == "Nelder-Mead")
# warning("one-diml optimization by Nelder-Mead is unreliable: use optimize")
#
- if(!HaveDriftHess){
+ if(!HaveDriftHess & (length(drift.par)>0)){
hess2 <- .Internal(optimhess(coef[drift.par], fDrift, NULL, conDrift))
HESS[drift.par,drift.par] <- hess2
}
- if(!HaveDiffHess){
+ if(!HaveDiffHess & (length(diff.par)>0)){
hess1 <- .Internal(optimhess(coef[diff.par], fDiff, NULL, conDiff))
HESS[diff.par,diff.par] <- hess1
}
@@ -351,14 +381,12 @@
solve(oout$hessian)
else matrix(numeric(0L), 0L, 0L)
-
-
- min <- oout$value
-
mycoef <- as.list(coef)
names(mycoef) <- nm
mycoef[fixed.par] <- fixed
+ min <- minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
+
new("mle", call = call, coef = coef, fullcoef = unlist(mycoef),
vcov = vcov, min = min, details = oout, minuslogl = minusquasilogl,
method = method)
@@ -390,9 +418,21 @@
diff.par <- yuima at model@parameter at diffusion
drift.par <- yuima at model@parameter at drift
- fullcoef <- c(diff.par, drift.par)
+
+ fullcoef <- NULL
+
+ if(length(diff.par)>0)
+ fullcoef <- diff.par
+
+ if(length(drift.par)>0)
+ fullcoef <- c(fullcoef, drift.par)
+
+# fullcoef <- c(diff.par, drift.par)
npar <- length(fullcoef)
-
+# cat("\nparam\n")
+# print(param)
+# cat("\nfullcoef\n")
+# print(fullcoef)
nm <- names(param)
oo <- match(nm, fullcoef)
if(any(is.na(oo)))
Added: pkg/yuima/man/CPoint.Rd
===================================================================
--- pkg/yuima/man/CPoint.Rd (rev 0)
+++ pkg/yuima/man/CPoint.Rd 2010-07-20 11:11:25 UTC (rev 109)
@@ -0,0 +1,115 @@
+\name{CPoint}
+\alias{CPoint}
+\alias{qmleL}
+\alias{qmleR}
+\title{Adaptive LASSO estimation for stochastic differential equations}
+\description{Adaptive LASSO estimation for stochastic differential equations.}
+\usage{
+CPoint(yuima, param1, param2, print=FALSE, plot=FALSE)
+qmleL(yuima, start, method="BFGS", t, print=FALSE, lower, upper, ...)
+qmleR(yuima, start, method="BFGS", t, print=FALSE, lower, upper, ...)
+}
+\arguments{
+ \item{yuima}{a yuima object.}
+ \item{param1}{parameter values before the change point t}
+ \item{param2}{parameter values after the change point t}
+ \item{print}{produces some output.}
+ \item{t}{right end point for \code{qmleL}, left end point for \code{qmleR}. See Details.}
+ \item{plot}{plot change point test statistics? Default is \code{FALSE}.}
+ \item{...}{passed to \code{\link{optim}} method. See Examples.}
+ \item{lower}{a named list for specifying lower bounds of parameters}
+ \item{upper}{a named list for specifying upper bounds of parameters}
+ \item{start}{initial values to be passed to the optimizer.}
+ \item{method}{see Details.}
+}
+\details{
+
+ \code{CPoint} estimates the change point using quasi-maximum likelihood approach.
+
+ Function \code{qmleL} estimates the parameters in the diffusion matrix using
+ observations up to time \code{t}.
+
+ Function \code{qmleR} estimates the parameters in the diffusion matrix using
+ observations from time \code{t} to the end.
+
+ Arguments in both \code{qmleL} and \code{qmleR} follow the same rules
+ as in \code{\link{qmle}}.
+
+
+}
+\value{
+ \item{ans}{a list with change point instant, and paramters before and after
+ the change point.}
+}
+\author{The YUIMA Project Team}
+\examples{
+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 at data@zoo.data[[1]]
+x1 <- as.numeric(x1[length(x1)])
+x2 <- yuima1 at 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 at data@zoo.data[[1]] <- c(yuima1 at data@zoo.data[[1]], yuima2 at data@zoo.data[[1]][-1])
+yuima at data@zoo.data[[2]] <- c(yuima1 at data@zoo.data[[2]], yuima2 at 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,low=low,method="L-BFGS-B")
+tmp1
+tmp2 <- qmleR(yuima,start=list(theta1.1=0.3,theta1.2=0.5), t=8.5,low=low,method="L-BFGS-B")
+tmp2
+
+
+# first stage changepoint estimator
+t.est2 <- CPoint(yuima,param1=tmp1$par,param2=tmp2$par)
+t.est2$tau
+
+
+# second stage estimation of parameters given first stage
+# change point estimator
+tmp11 <- qmleL(yuima,start=as.list(tmp1$par), t=t.est2$tau-0.1,method="L-BFGS-B")
+tmp11
+
+tmp21 <- qmleR(yuima,start=as.list(tmp2$par), t=t.est2$tau+0.1,method="L-BFGS-B")
+tmp21
+
+# second stage estimator of the change point
+CPoint(yuima,param1=tmp11$par,param2=tmp21$par)
+}
+
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
+\keyword{ts}
More information about the Yuima-commits
mailing list