[Yuima-commits] r265 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 26 17:07:31 CET 2013
Author: lorenzo
Date: 2013-12-26 17:07:31 +0100 (Thu, 26 Dec 2013)
New Revision: 265
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/R/qmle.R
pkg/yuima/R/toLatex.R
pkg/yuima/man/qmle.Rd
Log:
Modified qmle for carma model driven by a brownian motion
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2013-11-27 10:29:43 UTC (rev 264)
+++ pkg/yuima/DESCRIPTION 2013-12-26 16:07:31 UTC (rev 265)
@@ -1,9 +1,9 @@
Package: yuima
Type: Package
Title: The YUIMA Project package (unstable version)
-Version: 0.1.218
-Date: 2013-11-27
-Depends: methods, zoo, stats4, utils
+Version: 0.1.219
+Date: 2013-12-26
+Depends: methods, zoo, stats4, utils, Matrix
Suggests: cubature, mvtnorm
Author: YUIMA Project Team.
Maintainer: Stefano M. Iacus <stefano.iacus at R-project.org>
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2013-11-27 10:29:43 UTC (rev 264)
+++ pkg/yuima/NAMESPACE 2013-12-26 16:07:31 UTC (rev 265)
@@ -3,9 +3,11 @@
##importFrom("stats", "end", "time", "start")
importFrom("graphics", "plot")
importFrom("zoo")
+importFrom("Matrix")
importFrom(stats, confint)
importFrom(stats4)
+
importFrom(utils, toLatex)
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2013-11-27 10:29:43 UTC (rev 264)
+++ pkg/yuima/R/qmle.R 2013-12-26 16:07:31 UTC (rev 265)
@@ -65,14 +65,6 @@
}
-
-
-
-
-
-
-
-
### I have rewritten qmle as a version of ml.ql
### This function has an interface more similar to mle.
### ml.ql is limited in that it uses fixed names for drift and diffusion
@@ -84,6 +76,9 @@
qmle <- function(yuima, start, method="BFGS", fixed = list(), print=FALSE,
lower, upper, joint=FALSE, ...){
+ if(is(yuima at model, "yuima.carma")&& length(yuima at model@info at scale.par)!=0){
+ method<-"L-BFGS-B"
+ }
call <- match.call()
if( missing(yuima))
@@ -96,8 +91,39 @@
if( missing(start) )
yuima.stop("Starting values for the parameters are missing.")
+ #14/12/2013 We modify the QMLE function when the model is a Carma(p,q).
+ # In this case we use a two step procedure:
+ # First) The Coefficient are obtained by QMLE computed using the Kalman Filter.
+ # Second) Using the result in Brockwell, Davis and Yang (2007) we retrieve
+ # the underlying Levy. The estimated increments are used to find the Lévy parameters.
+
+# if(is(yuima at model, "yuima.carma")){
+# yuima.warm("two step procedure for carma(p,q)")
+# return(null)
+# }
+#
+
diff.par <- yuima at model@parameter at diffusion
- drift.par <- yuima at model@parameter at drift
+
+# 24/12
+ if(is(yuima at model, "yuima.carma") && length(diff.par)==0
+ && length(yuima at model@parameter at jump)!=0){
+ diff.par<-yuima at model@parameter at jump
+ }
+
+ if(is(yuima at model, "yuima.carma") && length(yuima at model@parameter at jump)!=0){
+ yuima.warn("carma(p,q): the qmle for a carma(p,q) driven by a Jump process will be implemented as soon as possible ")
+ return(null)
+ }
+
+ # 24/12
+ if(is(yuima at model, "yuima.carma") && length(yuima at model@info at lin.par)>0){
+ yuima.warn("carma(p,q): the case of lin.par will be implemented as soon as")
+ return(null)
+ }
+
+
+ drift.par <- yuima at model@parameter at drift
jump.par <- yuima at model@parameter at jump
measure.par <- yuima at model@parameter at measure
common.par <- yuima at model@parameter at common
@@ -107,12 +133,20 @@
JointOptim <- TRUE
yuima.warn("Drift and diffusion parameters must be different. Doing
joint estimation, asymptotic theory may not hold true.")
+ # 24/12
+# if(is(yuima at model, "yuima.carma")){
+# JointOptim <- TRUE
+# yuima.warm("Carma(p.q): The case of common parameters in Drift and Diffusion Term will be implemented as soon as possible,")
+# #return(null)
+# }
}
-
- if(length(jump.par)+length(measure.par)>0)
- yuima.stop("Cannot estimate the jump models, yet")
-
+ if(!is(yuima at model, "yuima.carma")){
+ if(length(jump.par)+length(measure.par)>0)
+ yuima.stop("Cannot estimate the jump models, yet")
+ }
+
+
if(!is.list(start))
yuima.stop("Argument 'start' must be of list type.")
@@ -166,13 +200,24 @@
d.size <- yuima at model@equation.number
+ if (is(yuima at model, "yuima.carma")){
+ # 24/12
+ d.size <-1
+ }
n <- length(yuima)[1]
env <- new.env()
assign("X", as.matrix(onezoo(yuima)), envir=env)
assign("deltaX", matrix(0, n-1, d.size), envir=env)
- assign("time", as.numeric(index(yuima at data@zoo.data[[1]])), envir=env)
- for(t in 1:(n-1))
+ if (is(yuima at model, "yuima.carma")){
+ #24/12 If we consider a carma model,
+ # the observations are only the first column of env$X
+ env$X<-as.matrix(env$X[,1])
+ env$deltaX<-as.matrix(env$deltaX[,1])
+ }
+ assign("time", as.numeric(index(yuima at data@zoo.data[[1]])), envir=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]]), envir=env)
@@ -197,7 +242,7 @@
}
oout <- NULL
- HESS <- matrix(0, length(nm), length(nm))
+ HESS <- matrix(0, length(nm), length(nm))
colnames(HESS) <- nm
rownames(HESS) <- nm
HaveDriftHess <- FALSE
@@ -207,6 +252,12 @@
if(length(start)>1){ #Â?multidimensional optim
oout <- optim(start, fj, method = method, hessian = TRUE, lower=lower, upper=upper)
HESS <- oout$hessian
+ if(is(yuima at model,"yuima.carma") && length(yuima at model@info at scale.par)!=0){
+ b0<-paste0(yuima at model@info at ma.par,"0",collapse="")
+ idx.b0<-match(b0,rownames(HESS))
+ HESS<-HESS[-idx.b0,]
+ HESS<-HESS[,-idx.b0]
+ }
HaveDriftHess <- TRUE
HaveDiffHess <- TRUE
} else { ### one dimensional optim
@@ -283,9 +334,15 @@
mydots$hessian <- FALSE
mydots$upper <- unlist( upper[ nm[idx.drift] ])
mydots$lower <- unlist( lower[ nm[idx.drift] ])
-
if(length(mydots$par)>1){
+ if(is(yuima at model, "yuima.carma")&& length(yuima at model@info at scale.par)!=0){
+ name_b0<-paste0(yuima at model@info at ma.par,"0",collapse="")
+ index_b0<-match(name_b0,nm)
+ mydots$lower[index_b0]<-1
+ mydots$upper[index_b0]<-1+10^(-7)
+ }
oout1 <- do.call(optim, args=mydots)
+ # oout1 <- optim(mydots$par,f,method = "L-BFGS-B" , lower = mydots$lower, upper = mydots$upper)
} else {
mydots$f <- mydots$fn
mydots$fn <- NULL
@@ -327,9 +384,10 @@
coef <- oout$par
control=list()
par <- coef
- names(par) <- c(diff.par, drift.par)
- nm <- c(diff.par, drift.par)
-
+ if(!is(yuima at model,"yuima.carma")){
+ names(par) <- c(diff.par, drift.par)
+ nm <- c(diff.par, drift.par)
+ }
# print(par)
# print(coef)
conDrift <- list(trace = 5, fnscale = 1,
@@ -370,7 +428,13 @@
if(!HaveDriftHess & (length(drift.par)>0)){
#hess2 <- .Internal(optimhess(coef[drift.par], fDrift, NULL, conDrift))
hess2 <- optimHess(coef[drift.par], fDrift, NULL, control=conDrift)
- HESS[drift.par,drift.par] <- hess2
+ HESS[drift.par,drift.par] <- hess2
+ if(is(yuima at model,"yuima.carma") && length(yuima at model@info at scale.par)!=0){
+ b0<-paste0(yuima at model@info at ma.par,"0",collapse="")
+ idx.b0<-match(b0,rownames(HESS))
+ HESS<-HESS[-idx.b0,]
+ HESS<-HESS[,-idx.b0]
+ }
}
if(!HaveDiffHess & (length(diff.par)>0)){
@@ -390,7 +454,13 @@
mycoef[fixed.par] <- fixed
min <- minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
-
+
+ dummycov<-matrix(0,length(coef),length(coef))
+ rownames(dummycov)<-names(coef)
+ colnames(dummycov)<-names(coef)
+ dummycov[rownames(vcov),colnames(vcov)]<-vcov
+ vcov<-dummycov
+
new("mle", call = call, coef = coef, fullcoef = unlist(mycoef),
vcov = vcov, min = min, details = oout, minuslogl = minusquasilogl,
method = method)
@@ -415,15 +485,24 @@
}
-
-
-
-
minusquasilogl <- function(yuima, param, print=FALSE, env){
diff.par <- yuima at model@parameter at diffusion
- drift.par <- yuima at model@parameter at drift
+ # 24/12
+ if(is(yuima at model, "yuima.carma") && length(yuima at model@info at lin.par)==0
+ && length(yuima at model@parameter at jump)!=0){
+ diff.par<-yuima at model@parameter at jump
+ }
+ # 24/12
+ if(is(yuima at model, "yuima.carma") && length(yuima at model@info at lin.par)>0 ){
+ yuima.warn("carma(p,q): the case of lin.par will be implemented as soon as")
+ return(null)
+ }
+
+
+ drift.par <- yuima at model@parameter at drift
+
fullcoef <- NULL
if(length(diff.par)>0)
@@ -457,49 +536,83 @@
n.theta <- n.theta1+n.theta2
d.size <- yuima at model@equation.number
+
+
n <- length(yuima)[1]
+
+ if (is(yuima at model, "yuima.carma")){
+ # 24/12
+ d.size <-1
+ # We build the two step procedure as described in
+ if(length(yuima at model@info at scale.par)!=0){
+ prova<-as.numeric(param)
+ names(prova)<-fullcoef[oo]
+ param<-prova[c(length(prova):1)]
+
+ y<-as.numeric(env$X)
+ tt<-env$time
+ p<-yuima at model@info at p
+ a<-param[c(1:p)]
+# a_names<-names(param[c(1:p)])
+# names(a)<-a_names
+ b<-param[c((p+1):(length(param)-p+1))]
+# b_names<-names(param[c((p+1):(length(param)-p+1))])
+# names(b)<-b_names
+ if(length(b)==1){
+ b<-1
+ } else{
+ indx_b0<-paste0(yuima at model@info at ma.par,"0",collapse="")
+ b[indx_b0]<-1
+ }
+ sigma<-tail(param,1)
+ strLog<-yuima.carma.loglik1(y, tt, a, b, sigma)
+ QL<-strLog$loglikCdiag
+ }else {
+ yuima.warn("carma(p,q): the scale parameter is equal to 1. We will implemented as soon as possible")
+ return(NULL)
+ }
+ } else{
+ drift <- drift.term(yuima, param, env)
+ diff <- diffusion.term(yuima, param, env)
+
+ QL <- 0
+
+ pn <- 0
+
+
+ vec <- env$deltaX-h*drift[-n,]
+
+ K <- -0.5*d.size * log( (2*pi*h) )
+
+ 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 <- K - 0.5*logdet-0.5*vec[t, ]^2/(h*yB)
+ QL <- QL+pn
+
+ }
+ } 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 <- log(1)
+ yuima.warn("singular diffusion matrix")
+ return(1e10)
+ }else{
+ pn <- K - 0.5*logdet +
+ ((-1/(2*h))*t(vec[t, ])%*%solve(yB)%*%vec[t, ])
+ QL <- QL+pn
+ }
+ }
+ }
+ }
- drift <- drift.term(yuima, param, env)
- diff <- diffusion.term(yuima, param, env)
- QL <- 0
-
- pn <- 0
-
-
- vec <- env$deltaX-h*drift[-n,]
-
- K <- -0.5*d.size * log( (2*pi*h) )
-
- 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 <- K - 0.5*logdet-0.5*vec[t, ]^2/(h*yB)
- QL <- QL+pn
-
- }
- } 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 <- log(1)
- yuima.warn("singular diffusion matrix")
- return(1e10)
- }else{
- pn <- K - 0.5*logdet +
- ((-1/(2*h))*t(vec[t, ])%*%solve(yB)%*%vec[t, ])
- QL <- QL+pn
- }
- }
- }
-
-
-
if(!is.finite(QL)){
yuima.warn("quasi likelihood is too small to calculate.")
return(1e10)
@@ -515,7 +628,134 @@
+MatrixA<-function (a)
+{
+ #Build Matrix A in the state space representation of Carma(p,q)
+ #given the autoregressive coefficient
+ pp = length(a)
+ af = cbind(rep(0, pp - 1), diag(pp - 1))
+ af = rbind(af, -a[pp:1])
+ return(af)
+}
+
+yuima.Vinfinity<-function(elForVInf,v){
+ # We find the infinity stationary variance-covariance matrix
+ A<-elForVInf$A
+ sigma<-elForVInf$sigma
+ p<-dim(A)[1]
+ matrixV<-matrix(0,p,p)
+ matrixV[upper.tri(matrixV,diag=TRUE)]<-v
+ matrixV[lower.tri(matrixV)]<-matrixV[upper.tri(matrixV)]
+ l<-rbind(matrix(rep(0,p-1),p-1,1),1)
+
+ RigSid<-l%*%t(l)
+ Matrixobj<-A%*%matrixV+matrixV%*%t(A)+sigma^2*RigSid
+ obj<-sum(Matrixobj^2)
+ obj
+}
+
+
+carma.kalman<-function(y, tt, p, q, a,bvector, sigma){
+
+ V_inf0<-matrix(diag(rep(1,p)),p,p)
+
+ A<-MatrixA(a)
+ u<-diff(tt)[1]
+ expA<-as.matrix(expm(A*u))
+ v<-as.numeric(V_inf0[upper.tri(V_inf0,diag=TRUE)])
+ elForVInf<-list(A=A,sigma=sigma)
+
+ V_inf_vect<-nlm(yuima.Vinfinity, v, elForVInf = elForVInf)$estimate
+ V_inf<-matrix(0,p,p)
+ V_inf[upper.tri(V_inf,diag=TRUE)]<-V_inf_vect
+ V_inf[lower.tri(V_inf)]<-V_inf[upper.tri(V_inf)]
+
+ V_inf[abs(V_inf)<=1.e-06]=0
+
+
+ SIGMA_err<-V_inf-expA%*%V_inf%*%t(expA)
+
+ statevar0<-matrix(rep(0, p),p,1)
+ Qmatr<-SIGMA_err
+
+ # set
+ statevar<-statevar0
+
+ # SigMatr<-expA%*%V_inf%*%t(expA)+Qmatr
+
+ SigMatr<-Qmatr
+ #SigMatr<-V_inf
+
+ zc<-matrix(bvector,1,p)
+ loglstar = 0
+ loglstar1 = 0
+ for(t in 1:length(tt)){
+ # prediction
+ statevar<-expA%*%statevar
+ SigMatr<-expA%*%SigMatr%*%t(expA)+Qmatr
+ # forecast
+ Uobs<-y[t]-zc%*%statevar
+ sd_2<-zc%*%SigMatr%*%t(zc)
+
+ #correction
+ Kgain<-SigMatr%*%t(zc)%*%solve(sd_2)
+ # SigMatr<-SigMatr-Kgain%*%as.numeric(sd_2)%*%t(Kgain)
+
+ statevar<-statevar+Kgain%*%Uobs
+ SigMatr<-SigMatr-Kgain%*%zc%*%SigMatr
+ # construction of log-likelihood
+ # loglstar1<-loglstar1+log(dnorm(as.numeric(Uobs), mean = 0, sd = sqrt(as.numeric(sd_2))))
+ # sdsig<-sqrt(as.numeric(sd_2))
+ term_int<--0.5*(log((2*pi)^(length(Uobs))*det(sd_2))+t(Uobs)%*%solve(sd_2)%*%Uobs)
+ loglstar<-loglstar+term_int
+ }
+ return(list(loglstar=as.numeric(loglstar),s2hat=as.numeric(sd_2)))
+}
+
+
+
+yuima.carma.loglik1<-function (y, tt, a, b, sigma)
+{
+ #This code compute the LogLik using kalman filter
+
+ # if(a_0!=0){we need to correct the Y_t for the mean}
+ # if(sigma!=1){we need to write}
+ p <- as.integer(length(a))
+
+ bvector <- rep(0, p)
+ q <- length(b)
+ bvector <- c(b, rep(0, p - q))
+
+
+ sigma<-sigma
+ y<-y
+
+ xxalt<-carma.kalman(y, tt, p, q, a,bvector,sigma)
+ list(loglikCdiag = xxalt$loglstar,s2hat=xxalt$s2hat)
+}
+
+
+loglik5 <- function(param) {
+ a<-param[1:pp]
+ b <- 1
+ if(qq>0){
+ b<-param[(pp + 1):(pp + qq)]
+ }
+
+ sigma<-tail(param,n=1)
+ xx <- yuima.carma.loglik1(y, tt, a, b, sigma)
+
+ return(xx$loglikCdiag)
+}
+
+
+
+
+
+
+
+
# returns the vector of log-transitions instead of the final quasilog
quasiloglvec <- function(yuima, param, print=FALSE, env){
Modified: pkg/yuima/R/toLatex.R
===================================================================
--- pkg/yuima/R/toLatex.R 2013-11-27 10:29:43 UTC (rev 264)
+++ pkg/yuima/R/toLatex.R 2013-12-26 16:07:31 UTC (rev 265)
@@ -8,8 +8,21 @@
mod <- object
if (class(object) == "yuima")
mod <- object at model
- if(class(mod) =="yuima.carma" && length(mod at info@lin.par)==0 )
- { yuima.warn("")
+ #if(class(mod) =="yuima.carma" && length(mod at info@lin.par)==0 )
+ if(class(mod) =="yuima.carma" )
+ {
+# yuima.warn("")
+
+
+ mysymb <- c("*", "alpha", "beta", "gamma", "delta", "rho",
+ "theta","sigma","mu", "sqrt")
+ # myrepl <- c(" \\cdot ", "\\alpha ", "\\beta ", "\\gamma ",
+ # "\\delta ", "\\rho ", "\\theta ", "\\sqrt ")
+ myrepl <- c(" \\cdot ", "\\alpha ", "\\beta ", "\\gamma ",
+ "\\delta ", "\\rho ", "\\theta ","\\sigma","\\mu", "\\sqrt ")
+ ns <- length(mysymb)
+
+
n.eq <- mod at equation.number
info <- mod at info
noise.var<-"W"
@@ -23,7 +36,7 @@
}
if(!length(info at loc.par)==0 && length(info at scale.par)==0){
- main.con<-paste(info at scale.par,"* \\ ", info at ma.par)
+ main.con<-paste(info at loc.par,"+ \\ ", info at ma.par)
}
if(!length(info at loc.par)==0 && !length(info at scale.par)==0){
@@ -44,6 +57,11 @@
"+ e",sprintf("d%s", noise.var),"\\left(",
mod at time.variable, "\\right) \\\\ \n")
dr<- paste(dr, "\\end{array}\\right.")
+#11/12
+ for (i in 1:ns) {
+ dr <- gsub(mysymb[i], myrepl[i], dr, fixed = "TRUE")
+ }
+
body <- paste("%%% Copy and paste the following output in your LaTeX file")
body <- c(body, paste("$$"))
body <- c(body, dr)
@@ -68,6 +86,11 @@
X<-paste(info at Latent.var,"(",mod at time.variable,")",
"=\\left[",latent.lab,
"\\right]'")
+
+ for (i in 1:ns) {
+ X <- gsub(mysymb[i], myrepl[i], X, fixed = "TRUE")
+ }
+
body <- c(body, X)
body <- c(body, paste("$$"))
# Vector Moving Average Coefficient.
@@ -75,16 +98,20 @@
#b.nozeros <-c(0:info at q)
- ma.lab0<-paste(paste(info at ma.par,0:(info at q),sep="_"),collapse=", \\ ")
+ # ma.lab0<-paste(paste(info at ma.par,0:(info at q),sep="_"),collapse=", \\ ")
+ ma.lab0<-paste(info at ma.par,0:(info at q),sep="_")
-
- if(length(ma.lab0)==1){ma.lab1<-ma.lab0}
- if(length(ma.lab0)==2){
- ma.lab0[1]<-paste(ma.lab0[1],",\\ ",sep="")
- # ma.lab0[2]<-paste(ma.lab0[2],"(",mod at time.variable,")",sep="")
- ma.lab1<-ma.lab0
- }
- if(length(ma.lab0)>2){
+ #if(length(ma.lab0)==1){ma.lab1<-ma.lab0}
+ if(info at q>=0 && info at q<=1){
+ ma.lab1<-paste(ma.lab0,collapse=", \\ ")}
+ #if(length(ma.lab0)==2){
+# if(info at q==1){
+# ma.lab0[1]<-paste(ma.lab0[1],",\\ ",sep="")
+# # ma.lab0[2]<-paste(ma.lab0[2],"(",mod at time.variable,")",sep="")
+# ma.lab1<-ma.lab0
+# }
+ #if(length(ma.lab0)>2){
+ if(info at q>1){
ma.lab1<-paste(ma.lab0[1],
",\\ ","\\ldots",
" \\ , \\ ",tail(ma.lab0,n=1))
@@ -102,6 +129,11 @@
ma.lab <- paste(ma.lab1," ,\\ 0, \\ \\ldots \\ , \\ 0")
}
Vector.ma <- paste(info at ma.par,"=","\\left[",ma.lab,"\\right]'")
+
+ for (i in 1:ns) {
+ Vector.ma <- gsub(mysymb[i], myrepl[i], Vector.ma, fixed = "TRUE")
+ }
+
body <- c(body, Vector.ma)
body <- c(body, paste("$$"))
@@ -110,9 +142,32 @@
noise.coef<-mod at diffusion
vect.e0 <- substr(tail(noise.coef,n=1), 13, nchar(tail(noise.coef,n=1)) -2)
+ vect.e1 <- substr(vect.e0, 2, nchar(vect.e0) -1)
+ dummy.e0<-strsplit(vect.e1,split="+",fixed=TRUE)
+
+
if (!length(mod at jump.variable)==0){
noise.coef <- mod at jump.coeff
vect.e0 <- substr(tail(noise.coef,n=1), 2, nchar(tail(noise.coef,n=1)) -1)
+ } else{
+ if(length(info at lin.par) != 0){
+
+ if (info at lin.par != info at ma.par){
+ dummy.e0<-c(dummy.e0[[1]][1], paste(dummy.e0[[1]][(2:length(dummy.e0[[1]]))],
+ "(",mod at time.variable,")",sep=""))
+ dummy.e0<-gsub(info at Latent.var,paste0(info at Latent.var,"_",collapse=""),dummy.e0)
+ dummy.e0<-gsub(info at lin.par,paste0(info at lin.par,"_",collapse=""),dummy.e0)
+ if(length(dummy.e0)>3){
+ vect.e0<-paste(dummy.e0[1],"+",dummy.e0[2],
+ "+ \\ldots +",tail(dummy.e0,n=1))
+ vect.e0<-paste("(",vect.e0,")")
+ } else{vect.e0<-paste("(",paste(dummy.e0,collapse="+"),")")}
+ }
+ # else{
+ # yuima.warm("The case of loc.par and scale.par will be implemented as soon as possible")
+ # return(NULL)
+ # }
+ }
}
if (info at p==1){vect.e <- vect.e0}
@@ -122,6 +177,13 @@
coeff.e<- paste("e","=","\\left[", vect.e , "\\right]'")
+ for (i in 1:ns) {
+ coeff.e <- gsub(mysymb[i], myrepl[i], coeff.e, fixed = "TRUE")
+ }
+
+
+
+
body <- c(body, coeff.e)
body <- c(body, paste("$$"))
# Matrix A
@@ -156,17 +218,18 @@
array.start<-paste0("\\begin{array}{",cent.col,"}\n",collapse="")
MATR.A<-paste("A ","=","\\left[",array.start, matrix.A, "\\end{array}\\right]'" )
+
+ for (i in 1:ns) {
+ MATR.A <- gsub(mysymb[i], myrepl[i], MATR.A, fixed = "TRUE")
+ }
+
+
body <- c(body, MATR.A)
body <- c(body, paste("$$"))
body <- structure(body, class = "Latex")
return(body)
- mysymb <- c("*", "alpha", "beta", "gamma", "delta", "rho",
- "theta","sigma","mu", "sqrt")
- # myrepl <- c(" \\cdot ", "\\alpha ", "\\beta ", "\\gamma ",
- # "\\delta ", "\\rho ", "\\theta ", "\\sqrt ")
- myrepl <- c(" \\cdot ", "\\alpha ", "\\beta ", "\\gamma ",
- "\\delta ", "\\rho ", "\\theta ","\\sigma","\\mu", "\\sqrt ")
+
} else{
n.eq <- mod at equation.number
dr <- paste("\\left(\\begin{array}{c}\n")
Modified: pkg/yuima/man/qmle.Rd
===================================================================
--- pkg/yuima/man/qmle.Rd 2013-11-27 10:29:43 UTC (rev 264)
+++ pkg/yuima/man/qmle.Rd 2013-12-26 16:07:31 UTC (rev 265)
@@ -1,198 +1,254 @@
-\name{qmle}
-%\alias{ql}
-%\alias{ml.ql}
-\alias{qmle}
-\alias{quasilogl}
-%\alias{LSE}
-%\alias{ml.ql2}
-\alias{rql}
-\alias{lse}
-%\alias{ql,ANY-method}
-%\alias{ml.ql,ANY-method}
-%\alias{ml.ql2,ANY-method}
-%\alias{rql,ANY-method}
-\title{Calculate quasi-likelihood and ML estimator of least squares estimator}
-\description{Calculate the quasi-likelihood and estimate of the parameters of the
- stochastic differential equation by the maximum likelihood method or least squares estimator
- of the drift parameter.}
-\usage{
-%ml.ql(yuima,theta2,theta1,h,theta2.lim=matrix(c(0,1),1,2),theta1.lim=matrix(c(0,1),1,2),print=FALSE,method,param,interval)
-%ql(yuima,theta2,theta1,h,print=FALSE,param)
-%rql(yuima,theta2,theta1,ptheta2,ptheta1,h,print=FALSE,param,prevparam)
-qmle(yuima, start, method="BFGS", fixed = list(), print=FALSE, lower, upper, joint=FALSE, ...)
-quasilogl(yuima, param, print=FALSE)
-}
-\arguments{
- \item{yuima}{a yuima object.}
-% \item{theta2,theta1}{parameters of the sdeModel.}
-% \item{h}{ time span of observations.}
-% \item{theta2.lim, theta1.lim}{matrixes to specify the domains of the
-% parameters. Vector can be available only if theta is a scalar.}
-% \item{ptheta2,ptheta1}{}
- \item{print}{you can see a progress of the estimation when print is TRUE.}
- \item{method}{see Details.}
- \item{param}{\code{list} of parameters for the quasi loglikelihood.}
-% \item{interval}{}
-% \item{prevparam}{}
- \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{fixed}{for conditional (quasi)maximum likelihood estimation.}
- \item{joint}{perform joint estimation or two stage estimation? by default \code{joint=FALSE}.}
- \item{...}{passed to \code{\link{optim}} method. See Examples.}
-}
-\details{
-% A function ql calculate the quasi-likelihood of a time series data X with any
-% parameters. A function ml.pl estimates parameters of the sdeModel by
-% maximizing the quasi-likelihood.
- \code{qmle} behaves more likely the standard \code{mle} function in \pkg{stats4} and
- argument \code{method} is one of the methods available in \code{\link{optim}}.
-
- \code{lse} calculates least squares estimators of the drift parameters. This is
- useful for initial guess of \code{qmle} estimation.
-
- \code{quasilogl} returns the valueof the quasi loglikelihood for a given
- \code{yuima} object and list of parameters \code{coef}.
-
-}
-\value{
- \item{QL}{a real value.}
- \item{opt}{a list with components the same as 'optim' function.}
-}
-\author{The YUIMA Project Team}
-\note{
- %The function ml.ql uses the function optim internally.
- The function qmle uses the function optim internally.
-}
-\examples{
-#dXt^e = -theta2 * Xt^e * dt + theta1 * dWt
-diff.matrix <- matrix(c("theta1"), 1, 1)
-ymodel <- setModel(drift=c("(-1)*theta2*x"), diffusion=diff.matrix, time.variable="t", state.variable="x", solve.variable="x")
-n <- 100
-
-ysamp <- setSampling(Terminal=(n)^(1/3), n=n)
-yuima <- setYuima(model=ymodel, sampling=ysamp)
-set.seed(123)
-yuima <- simulate(yuima, xinit=1, true.parameter=list(theta1=0.3,
-theta2=0.1))
-QL <- quasilogl(yuima, param=list(theta2=0.8, theta1=0.7))
-##QL <- ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)))
-QL
-
-## another way of parameter specification
-##param <- list(theta2=0.8, theta1=0.7)
-##QL <- ql(yuima, h=1/((n)^(2/3)), param=param)
-##QL
-
-
-## old code
-##system.time(
-##opt <- ml.ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)), c(0, 1), c(0, 1))
-##)
-##cat(sprintf("\nTrue param. theta2 = .3, theta1 = .1\n"))
-##print(coef(opt))
-
-
-system.time(
-opt2 <- qmle(yuima, start=list(theta1=0.8, theta2=0.7), lower=list(theta1=0,theta2=0),
- upper=list(theta1=1,theta2=1), method="L-BFGS-B")
-)
-cat(sprintf("\nTrue param. theta2 = .3, theta1 = .1\n"))
-print(coef(opt2))
-
-## initial guess for theta2 by least squares estimator
-tmp <- lse(yuima, start=list(theta2=0.7), lower=list(theta2=0), upper=list(theta2=1))
-tmp
-
-system.time(
-opt3 <- qmle(yuima, start=list(theta1=0.8, theta2=tmp), lower=list(theta1=0,theta2=0),
- upper=list(theta1=1,theta2=1), method="L-BFGS-B")
-)
-cat(sprintf("\nTrue param. theta2 = .3, theta1 = .1\n"))
-print(coef(opt3))
-
-
-## perform joint estimation? Non-optimal, just for didactic purposes
-system.time(
-opt4 <- qmle(yuima, start=list(theta1=0.8, theta2=0.7), lower=list(theta1=0,theta2=0),
- upper=list(theta1=1,theta2=1), method="L-BFGS-B", joint=TRUE)
-)
-cat(sprintf("\nTrue param. theta2 = .3, theta1 = .1\n"))
-print(coef(opt4))
-
-## old code
-##system.time(
-##opt <- ml.ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)), c(0, 1), c(0, 1), method="Newton")
-##)
-##cat(sprintf("\nTrue param. theta2 = .3, theta1 = .1\n"))
-##print(coef(opt))
-
-
-
-
-###multidimension case
-##dXt^e = - drift.matrix * Xt^e * dt + diff.matrix * dWt
-diff.matrix <- matrix(c("theta1.1","theta1.2", "1", "1"), 2, 2)
-
-drift.c <- c("-theta2.1*x1", "-theta2.2*x2", "-theta2.2", "-theta2.1")
-drift.matrix <- matrix(drift.c, 2, 2)
-
-ymodel <- setModel(drift=drift.matrix, diffusion=diff.matrix, time.variable="t",
- state.variable=c("x1", "x2"), solve.variable=c("x1", "x2"))
-n <- 100
-ysamp <- setSampling(Terminal=(n)^(1/3), n=n)
-yuima <- setYuima(model=ymodel, sampling=ysamp)
-set.seed(123)
-
-##xinit=c(x1, x2) #true.parameter=c(theta2.1, theta2.2, theta1.1, theta1.2)
-yuima <- simulate(yuima, xinit=c(1, 1), true.parameter=list(theta2.1=0.5, theta2.2=0.3, theta1.1=0.6, theta1.2=0.2))
-
-## theta2 <- c(0.8, 0.2) #c(theta2.1, theta2.2)
-##theta1 <- c(0.7, 0.1) #c(theta1.1, theta1.2)
-## QL <- ql(yuima, theta2, theta1, h=1/((n)^(2/3)))
-## QL
-
-## ## another way of parameter specification
-## #param <- list(theta2=theta2, theta1=theta1)
-## #QL <- ql(yuima, h=1/((n)^(2/3)), param=param)
-## #QL
-
-## theta2.1.lim <- c(0, 1)
-## theta2.2.lim <- c(0, 1)
-## theta1.1.lim <- c(0, 1)
-## theta1.2.lim <- c(0, 1)
-## theta2.lim <- t( matrix( c(theta2.1.lim, theta2.2.lim), 2, 2) )
-## theta1.lim <- t( matrix( c(theta1.1.lim, theta1.2.lim), 2, 2) )
-
-## system.time(
-## opt <- ml.ql(yuima, theta2, theta1, h=1/((n)^(2/3)), theta2.lim, theta1.lim)
-## )
-## opt at coef
-
-system.time(
-opt2 <- qmle(yuima, start=list(theta2.1=0.8, theta2.2=0.2, theta1.1=0.7, theta1.2=0.1),
- lower=list(theta1.1=.1,theta1.2=.1,theta2.1=.1,theta2.2=.1),
- upper=list(theta1.1=4,theta1.2=4,theta2.1=4,theta2.2=4), method="L-BFGS-B")
-)
-opt2 at coef
-summary(opt2)
-
-## unconstrained optimization
-system.time(
-opt3 <- qmle(yuima, start=list(theta2.1=0.8, theta2.2=0.2, theta1.1=0.7, theta1.2=0.1))
-)
-opt3 at coef
-summary(opt3)
-
-
-quasilogl(yuima, param=list(theta2.1=0.8, theta2.2=0.2, theta1.1=0.7, theta1.2=0.1))
-
-##system.time(
-##opt <- ml.ql(yuima, theta2, theta1, h=1/((n)^(2/3)), theta2.lim, theta1.lim, method="Newton")
-##)
-##opt at coef
-##}
-
-% Add one or more standard keywords, see file 'KEYWORDS' in the
-% R documentation directory.
-\keyword{ts}
+\name{qmle}
+%\alias{ql}
+%\alias{ml.ql}
+\alias{qmle}
+\alias{quasilogl}
+%\alias{LSE}
+%\alias{ml.ql2}
+\alias{rql}
+\alias{lse}
+%\alias{ql,ANY-method}
+%\alias{ml.ql,ANY-method}
+%\alias{ml.ql2,ANY-method}
+%\alias{rql,ANY-method}
+\title{Calculate quasi-likelihood and ML estimator of least squares estimator}
+\description{Calculate the quasi-likelihood and estimate of the parameters of the
+ stochastic differential equation by the maximum likelihood method or least squares estimator
+ of the drift parameter.}
+\usage{
+%ml.ql(yuima,theta2,theta1,h,theta2.lim=matrix(c(0,1),1,2),theta1.lim=matrix(c(0,1),1,2),print=FALSE,method,param,interval)
+%ql(yuima,theta2,theta1,h,print=FALSE,param)
+%rql(yuima,theta2,theta1,ptheta2,ptheta1,h,print=FALSE,param,prevparam)
+qmle(yuima, start, method="BFGS", fixed = list(), print=FALSE, lower, upper, joint=FALSE, ...)
+quasilogl(yuima, param, print=FALSE)
+}
+\arguments{
+ \item{yuima}{a yuima object.}
+% \item{theta2,theta1}{parameters of the sdeModel.}
+% \item{h}{ time span of observations.}
+% \item{theta2.lim, theta1.lim}{matrixes to specify the domains of the
+% parameters. Vector can be available only if theta is a scalar.}
+% \item{ptheta2,ptheta1}{}
+ \item{print}{you can see a progress of the estimation when print is TRUE.}
+ \item{method}{see Details.}
+ \item{param}{\code{list} of parameters for the quasi loglikelihood.}
+% \item{interval}{}
+% \item{prevparam}{}
+ \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{fixed}{for conditional (quasi)maximum likelihood estimation.}
+ \item{joint}{perform joint estimation or two stage estimation? by default \code{joint=FALSE}.}
+ \item{...}{passed to \code{\link{optim}} method. See Examples.}
+}
+\details{
+% A function ql calculate the quasi-likelihood of a time series data X with any
+% parameters. A function ml.pl estimates parameters of the sdeModel by
+% maximizing the quasi-likelihood.
+ \code{qmle} behaves more likely the standard \code{mle} function in \pkg{stats4} and
+ argument \code{method} is one of the methods available in \code{\link{optim}}.
+
+ \code{lse} calculates least squares estimators of the drift parameters. This is
+ useful for initial guess of \code{qmle} estimation.
+
+ \code{quasilogl} returns the valueof the quasi loglikelihood for a given
+ \code{yuima} object and list of parameters \code{coef}.
+
+}
+\value{
+ \item{QL}{a real value.}
+ \item{opt}{a list with components the same as 'optim' function.}
+}
+\author{The YUIMA Project Team}
+\note{
+ %The function ml.ql uses the function optim internally.
+ The function qmle uses the function optim internally.
+}
+\examples{
+#dXt^e = -theta2 * Xt^e * dt + theta1 * dWt
+diff.matrix <- matrix(c("theta1"), 1, 1)
+ymodel <- setModel(drift=c("(-1)*theta2*x"), diffusion=diff.matrix, time.variable="t", state.variable="x", solve.variable="x")
+n <- 100
+
+ysamp <- setSampling(Terminal=(n)^(1/3), n=n)
+yuima <- setYuima(model=ymodel, sampling=ysamp)
+set.seed(123)
+yuima <- simulate(yuima, xinit=1, true.parameter=list(theta1=0.3,
+theta2=0.1))
+QL <- quasilogl(yuima, param=list(theta2=0.8, theta1=0.7))
+##QL <- ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)))
+QL
+
+## another way of parameter specification
+##param <- list(theta2=0.8, theta1=0.7)
+##QL <- ql(yuima, h=1/((n)^(2/3)), param=param)
+##QL
+
+
+## old code
+##system.time(
+##opt <- ml.ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)), c(0, 1), c(0, 1))
+##)
+##cat(sprintf("\nTrue param. theta2 = .3, theta1 = .1\n"))
+##print(coef(opt))
+
+
+system.time(
+opt2 <- qmle(yuima, start=list(theta1=0.8, theta2=0.7), lower=list(theta1=0,theta2=0),
+ upper=list(theta1=1,theta2=1), method="L-BFGS-B")
+)
+cat(sprintf("\nTrue param. theta2 = .3, theta1 = .1\n"))
+print(coef(opt2))
+
+## initial guess for theta2 by least squares estimator
+tmp <- lse(yuima, start=list(theta2=0.7), lower=list(theta2=0), upper=list(theta2=1))
+tmp
+
+system.time(
+opt3 <- qmle(yuima, start=list(theta1=0.8, theta2=tmp), lower=list(theta1=0,theta2=0),
+ upper=list(theta1=1,theta2=1), method="L-BFGS-B")
+)
+cat(sprintf("\nTrue param. theta2 = .3, theta1 = .1\n"))
+print(coef(opt3))
+
+
+## perform joint estimation? Non-optimal, just for didactic purposes
+system.time(
+opt4 <- qmle(yuima, start=list(theta1=0.8, theta2=0.7), lower=list(theta1=0,theta2=0),
+ upper=list(theta1=1,theta2=1), method="L-BFGS-B", joint=TRUE)
+)
+cat(sprintf("\nTrue param. theta2 = .3, theta1 = .1\n"))
+print(coef(opt4))
+
+## old code
+##system.time(
+##opt <- ml.ql(yuima, 0.8, 0.7, h=1/((n)^(2/3)), c(0, 1), c(0, 1), method="Newton")
+##)
+##cat(sprintf("\nTrue param. theta2 = .3, theta1 = .1\n"))
+##print(coef(opt))
+
+
+
+
+###multidimension case
+##dXt^e = - drift.matrix * Xt^e * dt + diff.matrix * dWt
+diff.matrix <- matrix(c("theta1.1","theta1.2", "1", "1"), 2, 2)
+
+drift.c <- c("-theta2.1*x1", "-theta2.2*x2", "-theta2.2", "-theta2.1")
+drift.matrix <- matrix(drift.c, 2, 2)
+
+ymodel <- setModel(drift=drift.matrix, diffusion=diff.matrix, time.variable="t",
+ state.variable=c("x1", "x2"), solve.variable=c("x1", "x2"))
+n <- 100
+ysamp <- setSampling(Terminal=(n)^(1/3), n=n)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 265
More information about the Yuima-commits
mailing list