[Yuima-commits] r5 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Oct 27 09:50:46 CET 2009
Author: iacus
Date: 2009-10-27 09:50:46 +0100 (Tue, 27 Oct 2009)
New Revision: 5
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/R/adaBayes.R
pkg/yuima/R/quasi-likelihood.R
pkg/yuima/R/rng.R
pkg/yuima/R/simulate.R
pkg/yuima/R/yuima.model.R
pkg/yuima/man/adaBayes.Rd
pkg/yuima/man/quasi-likelihood.Rd
pkg/yuima/man/simulate.Rd
Log:
last update from csv
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2009-10-07 12:09:39 UTC (rev 4)
+++ pkg/yuima/DESCRIPTION 2009-10-27 08:50:46 UTC (rev 5)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package
-Version: 0.0.64
-Date: 2009-09-13
+Version: 0.0.66
+Date: 2009-10-25
Depends: methods, zoo
Author: Yuima project team.
Maintainer: Who to complain to <yourfault at somewhere.net>
Modified: pkg/yuima/R/adaBayes.R
===================================================================
--- pkg/yuima/R/adaBayes.R 2009-10-07 12:09:39 UTC (rev 4)
+++ pkg/yuima/R/adaBayes.R 2009-10-27 08:50:46 UTC (rev 5)
@@ -23,7 +23,7 @@
tmp <- ((n.size*d.size/2) * log( (2*pi*h) ) + ql(yuima,theta2=theta2,theta1=theta1,h=h))
return(tmp)
}
-
+
## rHn function
rHn <- function(yuima,h,theta1,theta2,theta2.offset,theta1.offset){
d.size <- yuima at model@equation.number
@@ -31,7 +31,7 @@
tmp <- ((n.size*d.size/2) * log( (2*pi*h) ) + rql(yuima,theta2=theta2,theta1=theta1,theta2.offset,theta1.offset,h=h))
return(tmp)
}
-
+
## Bayes estimator for theta1
bayes.theta1.tilda <- function(yuima,prior1=prior1){
rHn.tmp <- function(yuima,h,theta1,theta2,theta2.offset,theta1.offset){
@@ -70,11 +70,11 @@
}
return( numerator$value/denominator$value )
}
-
+
##Bayes estimator for theta2
bayes.theta2.tilda <- function(yuima,prior2=prior2){
Hn.tmp <- function(yuima,h,theta1,theta2){
- tmp <- numeric( length(theta2) )
+ tmp <- numeric(length(theta2))
for(i in 1:length(theta2)){
Hn1 <- Hn( yuima , h , theta1 , theta2[i] )
tmp[i] <- exp( Hn1 -Hn2 )
Modified: pkg/yuima/R/quasi-likelihood.R
===================================================================
--- pkg/yuima/R/quasi-likelihood.R 2009-10-07 12:09:39 UTC (rev 4)
+++ pkg/yuima/R/quasi-likelihood.R 2009-10-27 08:50:46 UTC (rev 5)
@@ -8,7 +8,7 @@
modelstate <- yuima at model@state.variable
modelpara <- yuima at model@parameter at diffusion
DIFFUSION <- yuima at model@diffusion
- division <- length(yuima)
+ division <- length(yuima)[1]
diff <- array(0,dim=c(d.size,r.size,division))
X <- as.matrix(onezoo(yuima))
for(i in 1:length(modelpara)){
@@ -36,7 +36,7 @@
modelstate <- yuima at model@state.variable
modelpara <- yuima at model@parameter at drift
DRIFT <- yuima at model@drift
- division <- length(yuima)
+ division <- length(yuima)[1]
drift <- matrix(0,division,d.size)
X <- as.matrix(onezoo(yuima))
for(i in 1:length(modelpara)){
@@ -99,7 +99,7 @@
d.size <- yuima at model@equation.number
- division <- length(yuima)
+ division <- length(yuima)[1]
X <- as.matrix(onezoo(yuima))
deltaX <- matrix(0,division-1,d.size)
for(t in 1:(division-1)){
@@ -118,12 +118,20 @@
if(det(as.matrix(B[,,t]))==0){
pn[t] <- log(1)
}else{
- pn[t] <- log(1/((2*pi*h)^(d.size/2)*det(as.matrix(B[,,t]))^(1/2)) * exp((-1/(2*h))*t(deltaX[t,]-h*drift[t,])%*%solve(as.matrix(B[,,t]))%*%(deltaX[t,]-h*drift[t,])))
+ pn[t] <- log(1/((2*pi*h)^(d.size/2)*det(as.matrix(B[,,t]))^(1/2)) *
+ exp((-1/(2*h))*t(deltaX[t,]-h*drift[t,])%*%solve(as.matrix(B[,,t]))%*%(deltaX[t,]-h*drift[t,])))
QL <- QL+pn[t]
+ if(pn[t]==-Inf && FALSE){
+ cat("t:",t, "\n")
+ cat("B[,,t]:",B[,,t], "\n")
+ cat("det(B):",det(as.matrix(B[,,t])),"\n")
+ cat("deltaX[t,]", deltaX[t,], "\n")
+ cat("drift[t,]", drift[t,], "\n")
+ }
}
}
if(QL==-Inf){
- warning("quasi likelihood is to small to calculate.")
+ warning("quasi likelihood is too small too calculate.")
}
if(print==TRUE){
print(paste("QL:",QL," theta2:",theta2," theta1:",theta1))
@@ -168,7 +176,7 @@
d.size <- yuima at model@equation.number
- division <- length(yuima)
+ division <- length(yuima)[1]
X <- as.matrix(onezoo(yuima))
deltaX <- matrix(0,division-1,d.size)
for(t in 1:(division-1)){
@@ -261,7 +269,7 @@
ql.opt <- function(theta=c(theta2,theta1)){
return(ql(yuima,theta2=theta[1:length(theta2)],theta1=theta[(length(theta2)+1):length(theta)],h=h,print=print))
}
-
+
opt <- constrOptim(c(theta2,theta1),ql.opt,NULL,ui=ui,ci=ci,control=list(fnscale=-1),outer.iterations=500)
if( opt$convergence != 0){
print("WARNING:optimization did not converge.")
Modified: pkg/yuima/R/rng.R
===================================================================
--- pkg/yuima/R/rng.R 2009-10-07 12:09:39 UTC (rev 4)
+++ pkg/yuima/R/rng.R 2009-10-27 08:50:46 UTC (rev 5)
@@ -6,7 +6,7 @@
if( gamma <= 0 )
stop("deltat must be positive value.")
V <- rchisq(x,df=1);
-
+
z1 <- ( delta/gamma + V/(2*gamma^2) ) - sqrt( V*delta/(gamma^3) + ( V/(2*gamma^2) )^2 )
U <- runif(x,min=0,max=1)
idx <- which( U < (delta/(delta+gamma*z1)) )
@@ -14,25 +14,65 @@
ret <- numeric(x)
ret[idx] <- z1[idx]
ret[-idx] <- z2
-
+
return(ret)
}
-rNIG <- function(x,alpha=1,beta=0,delta=1,mu=0){
- gamma <- sqrt(alpha^2 - beta^2)
- if (gamma == 0) {
- V = rnorm(x)^2
- Z = delta * delta/V
- X = sqrt(Z) * rnorm(x)
- }else{
- Z <- rIG(x,delta,gamma)
- N <- rnorm(x,0,1)
- X <- mu + beta*Z + sqrt(Z)*N
+## multi-variate normal inverse Gaussian
+rNIG <- function(x,alpha=1,beta=0,delta=1,mu=0,Lambda){
+
+ ## Error check
+ if(length(mu)!=length(beta)){
+ stop("Error: wrong input dimension.")
}
+ if(missing(Lambda)){
+ ## univariate case
+ gamma <- sqrt(as.numeric(alpha^2 - beta^2))
+ }else{
+ if( det(Lambda) < 10^(-15)){
+ stop("Determinant of Lambda must be one.")
+ }
+ gamma <- sqrt(as.numeric(alpha^2 - t(beta) %*% Lambda %*% beta))
+ }
+
+ if (gamma == 0)
+ stop("alpha^2 - beta^2*mu must be positive value.")
+
+ tau <- rIG(x,delta,gamma)
+ eta <- rnorm(length(beta))
+ if(missing(Lambda)){ ## uni-variate
+ X <- mu + beta*tau + sqrt(tau)*eta
+ }else{ ## multi-variate
+ sqrt.L <- svd(Lambda)
+ sqrt.L <- sqrt.L$u %*% diag(sqrt(sqrt.L$d)) %*% t(sqrt.L$v)
+
+ z <- mu + t(matrix(rep(tau,length(beta)),x,length(beta))) * matrix(rep(Lambda %*% beta,x),length(beta),x)
+ + t(matrix(rep(tau,length(beta)),x,length(beta))) * matrix(rep(Lambda %*% eta,x),length(beta),x)
+
+ X <- z
+ }
+
return(X)
}
+
+## rNIG <- function(x,alpha=1,beta=0,delta=1,mu=0){
+## gamma <- sqrt(alpha^2 - beta^2)
+## if (gamma == 0) {
+## V = rnorm(x)^2
+## Z = delta * delta/V
+## X = sqrt(Z) * rnorm(x)
+## }else{
+## Z <- rIG(x,delta,gamma)
+## N <- rnorm(x,0,1)
+## X <- mu + beta*Z + sqrt(Z)*N
+## }
+## return(X)
+## }
+
+
+
rbgamma <- function(x,delta.plus,gamma.plus,delta.minus,gamma.minus){
if( delta.plus <= 0 )
stop("delta.plus must be positive value.")
@@ -47,31 +87,92 @@
}
-## temporaly, desined only for univariate.
+## temporaly, desined only for univariate.
rngamma <- function(x,lambda,alpha,beta,mu,Lambda){
- tmp <- alpha^2 - t(beta) %*% Lambda %*% beta
+ ## Error check
+ if(length(mu)!=length(beta)){
+ stop("Error: wrong input dimension.")
+ }
+ if(missing(Lambda)){
+ ## univariate case
+ tmp <- as.numeric(alpha^2 - beta^2)
+ }else{
+ if( det(Lambda) < 10^(-15)){
+ stop("Determinant of Lambda must be one.")
+ }
+ tmp <- as.numeric(alpha^2 - t(beta) %*% Lambda %*% beta)
+ }
+
if( lambda <= 0 )
stop("lambda must be positive value.")
if( alpha <= 0 )
stop("alpha must be positive value.")
- if( tmp <= 0 )
+ if( tmp <=0)
stop("alpha^2 - beta^2*mu must be positive value.")
-# if( det(Lambda) != 1)
- if( Lambda !=1)
- stop("Determinant of Lambda must be one.")
-
+
tau <- rgamma(x,lambda,tmp/2)
- eta <- rnorm(x)
- z <- mu + beta * tau * Lambda + sqrt(tau * Lambda) * eta
- X <- z
+ eta <- rnorm(length(beta))
+ if(missing(Lambda)){
+ z <- mu + beta * tau + sqrt(tau) * eta
+ X <- z
+ }else{
+ sqrt.L <- svd(Lambda)
+ sqrt.L <- sqrt.L$u %*% diag(sqrt(sqrt.L$d)) %*% t(sqrt.L$v)
+
+ z <- mu + t(matrix(rep(tau,length(beta)),x,length(beta))) * matrix(rep(Lambda %*% beta,x),length(beta),x)
+ + t(matrix(rep(tau,length(beta)),x,length(beta))) * matrix(rep(Lambda %*% eta,x),length(beta),x)
+
+ X <- z
+ }
return(X)
}
-rstable <- function(x,alpha,beta,sigma, gamma, eps){
+## ## temporaly, desined only for univariate.
+## rngamma <- function(x,lambda,alpha,beta,mu){
+## tmp <- alpha^2 - beta^2
+## if( lambda <= 0 )
+## stop("lambda must be positive value.")
+## if( alpha <= 0 )
+## stop("alpha must be positive value.")
+## if( tmp <= 0 )
+## stop("alpha^2 - beta^2*mu must be positive value.")
+## # if( det(Lambda) != 1)
+## ## if( Lambda !=1)
+## ## stop("Determinant of Lambda must be one.")
+
+## tau <- rgamma(x,lambda,tmp/2)
+## eta <- rnorm(x)
+## ## z <- mu + beta * tau * Lambda + sqrt(tau * Lambda) * eta
+## z <- mu + beta * tau + sqrt(tau) * eta
+## X <- z
+## return(X)
+## }
+
+## rngamma <- function(x,lambda,alpha,beta,mu,Lambda){
+## tmp <- alpha^2 - t(beta) %*% Lambda %*% beta
+## if( lambda <= 0 )
+## stop("lambda must be positive value.")
+## if( alpha <= 0 )
+## stop("alpha must be positive value.")
+## if( tmp <= 0 )
+## stop("alpha^2 - beta^2*mu must be positive value.")
+## # if( det(Lambda) != 1)
+## if( Lambda !=1)
+## stop("Determinant of Lambda must be one.")
+
+## tau <- rgamma(x,lambda,tmp/2)
+## eta <- rnorm(x)
+## z <- mu + beta * tau * Lambda + sqrt(tau * Lambda) * eta
+## X <- z
+## return(X)
+## }
+
+
+rstable <- function(x,alpha,beta,sigma, gamma){
a <- (1 + (beta*tan(alpha*pi/2))^2)^(1/(2*alpha))
b <- atan(beta*tan(alpha*pi/2))/alpha
-
+
u <- runif(x,-pi/2,pi/2)
v <- rexp(x,1)
@@ -81,6 +182,24 @@
s <- (2/pi) * ((pi/2 +beta*u)*tan(u) - beta * log(v*cos(u)/(beta*u + pi/2)))
}
- X <- (eps^(1/alpha)) * sigma * s + gamma * eps * rep(1,x)
+## X <- (eps^(1/alpha)) * sigma * s + gamma * eps * rep(1,x)
+ X <- sigma * s + gamma * rep(1,x)
return(X)
}
+
+## rstable <- function(x,alpha,beta,sigma, gamma, eps){
+## a <- (1 + (beta*tan(alpha*pi/2))^2)^(1/(2*alpha))
+## b <- atan(beta*tan(alpha*pi/2))/alpha
+
+## u <- runif(x,-pi/2,pi/2)
+## v <- rexp(x,1)
+
+## if(alpha!=1){
+## s <- a * (sin(alpha*(u+b))/cos(u)^(1/alpha)) * (cos(u-alpha*(u+b))/v)^((1-alpha)/alpha)
+## }else{
+## s <- (2/pi) * ((pi/2 +beta*u)*tan(u) - beta * log(v*cos(u)/(beta*u + pi/2)))
+## }
+
+## X <- (eps^(1/alpha)) * sigma * s + gamma * eps * rep(1,x)
+## return(X)
+## }
Modified: pkg/yuima/R/simulate.R
===================================================================
--- pkg/yuima/R/simulate.R 2009-10-07 12:09:39 UTC (rev 4)
+++ pkg/yuima/R/simulate.R 2009-10-27 08:50:46 UTC (rev 5)
@@ -22,10 +22,10 @@
##:: simulate
##:: solves SDE and returns result
setGeneric("simulate",
- function(yuima, xinit, true.parameter, space.discretized=FALSE, increment)
+ function(yuima, xinit, true.parameter, space.discretized=FALSE, increment.W=NULL, increment.L=NULL)
standardGeneric("simulate")
)
-setMethod("simulate", "yuima", function(yuima, xinit, true.parameter, space.discretized=FALSE, increment){
+setMethod("simulate", "yuima", function(yuima, xinit, true.parameter, space.discretized=FALSE, increment.W=NULL, increment.L=NULL){
##:: error checks
if(missing(yuima)){
cat("\nyuima object is missing.\n")
@@ -42,7 +42,7 @@
## readline()
if(missing(true.parameter)){
- true.parameter <- numeric(length(sdeModel at parameter@all))
+ true.parameter <- as.list(numeric(length(sdeModel at parameter@all)))
}
r.size <- sdeModel at noise.number
@@ -72,20 +72,34 @@
}
##:: Error check for increment specified version.
- if(!missing(increment)){
+ if(!missing(increment.W)){
if(space.discretized == TRUE){
cat("\nParameter increment must be invalid if space.discretized=TRUE.\n")
return(NULL)
- }else if(dim(increment)[1] != r.size){
+ }else if(dim(increment.W)[1] != r.size){
cat("\nLength of increment's row must be same as yuima at model@noise.number.\n")
return(NULL)
- }else if(dim(increment)[2] != division){
+ }else if(dim(increment.W)[2] != division){
cat("\nLength of increment's column must be same as yuima at sampling@division[1].\n")
return(NULL)
}
}
+ ##:: Error check for increment specified version.
+ if(!missing(increment.L)){
+ if(space.discretized == TRUE){
+ cat("\nParameter increment must be invalid if space.discretized=TRUE.\n")
+ return(NULL)
+ }else if(dim(increment.L)[1] != r.size){
+ cat("\nLength of increment's row must be same as yuima at model@noise.number.\n")
+ return(NULL)
+ }else if(dim(increment.L)[2] != division){
+ cat("\nLength of increment's column must be same as yuima at sampling@division[1].\n")
+ return(NULL)
+ }
+ }
+
##:: check if DRIFT and/or DIFFUSION has values
has.drift <- sum(as.character(sdeModel at drift) != "(0)")
var.in.diff <- is.logical(any(match(unlist(lapply(sdeModel at diffusion, all.vars)), sdeModel at state.variable)))
@@ -95,15 +109,15 @@
modeltime <- sdeModel at time.variable
V0 <- sdeModel at drift
V <- sdeModel at diffusion
-
+
par.len <- length(sdeModel at parameter@all)
if(par.len>0){
for(i in 1:par.len){
pars <- sdeModel at parameter@all[i]
- assign(pars, true.parameter[i])
+ assign(pars, true.parameter[[i]])
}
}
-
+
##:: Initialization
##:: set time step
delta <- Terminal/division
@@ -207,12 +221,15 @@
##:: using Euler-Maruyama method
##:: Diffusion terms
- if( missing(increment)){
+ if( missing(increment.W)){
dW <- rnorm(division * r.size, 0, sqrt(delta))
dW <- t(matrix(dW, nrow=division, ncol=r.size))
}else{
- dW <- increment
+ dW <- increment.W
}
+
+ ## [TBC] Levy increment¤¬»ØÄꤵ¤ì¤¿¾ì¹ç¤Ë¤â, ¥·¥ß¥å¥ì¡¼¥·¥ç¥ó¤òincrement¤òÍѤ¤¤Æ¹Ô¤¦¤è¤¦¤Ë.
+
if(has.drift){ ##:: consider drift term to be one of the diffusion term(dW=1)
dW <- rbind( rep(1, division)*delta , dW)
}
@@ -308,7 +325,7 @@
##:: make expression to create iid rand J
if(grep("^[dexp|dnorm|dgamma]", sdeModel at measure$df$expr)) {
##:: e.g. dnorm(z,1,1) -> rnorm(mu.size*N_sharp,1,1)
- F <- parse(text=gsub("^d(.+?)\\(.+?,", "r\\1(mu.size*N_sharp,", sdeModel at measure$df$expr, perl=TRUE))
+ F <- suppressWarnings(parse(text=gsub("^d(.+?)\\(.+?,", "r\\1(mu.size*N_sharp,", sdeModel at measure$df$expr, perl=TRUE)))
}else{
stop("Sorry. CP only supports dexp, dnorm and dgamma yet.")
}
@@ -333,23 +350,29 @@
}else if(sdeModel at measure.type=="code"){ ##:: code type
##:: Jump terms
- code <- sub("^(.+?)\\(.+", "\\1", sdeModel at measure$df$expr, perl=TRUE)
- args <- unlist(strsplit(sub("^.+?\\((.+)\\)", "\\1", sdeModel at measure$df$expr, perl=TRUE), ","))
-
+ code <- suppressWarnings(sub("^(.+?)\\(.+", "\\1", sdeModel at measure$df$expr, perl=TRUE))
+ args <- eval(parse(text=paste("list(",suppressWarnings(sub("^.+?\\(.+?,(.+)", "\\1", sdeModel at measure$df$expr, perl=TRUE)))))
+
+ # args <- unlist(strsplit(suppressWarnings(sub("^.+?\\((.+)\\)", "\\1", sdeModel at measure$df$expr, perl=TRUE)), ","))
+ #print(args);readline()
dZ <- switch(code,
- rNIG=paste("rNIG(division, ", args[2], ", ", args[3], ", ", args[4], "*delta, ", args[5], "*delta)"),
- rIG=paste("rIG(division,", args[2], "*delta, ", args[3], ")"),
- rgamma=paste("rgamma(division, ", args[2], "*delta, ", args[3], ")"),
- rbgamma=paste("rbgamma(division, ", args[2], "*delta, ", args[3], ", ", args[4], "*delta, ", args[5], ")"),
- rngamma=paste("rngamma(division, ", args[2], "*delta, ", args[3], ", ", args[4], ", ", args[5], "*delta, ", args[6], ")"),
- rstable=paste("rstable(division, ", args[2], ", ", args[3], ", ", args[4], ", ", args[5], ", ", args[6], ")")
+ #rNIG=paste("rNIG(division, ", args[2], ", ", args[3], ", ", args[4], "*delta, ", args[5], "*delta)"),
+ rNIG=rNIG(division, args[[1]], args[[2]], args[[3]]*delta, args[[4]]*delta,args[[5]]*delta),
+ rIG=rIG(division,args[[1]]*delta, args[[2]]),
+ rgamma=rgamma(division,args[[1]]*delta, args[[2]]),
+ rbgamma=rbgamma(division, args[[1]]*delta, args[[2]], args[[3]]*delta, args[[4]]),
+## rngamma=paste("rngamma(division, ", args[2], "*delta, ", args[3], ", ", args[4], ", ", args[5], "*delta, ", args[6], ")"),
+ rngamma=rngamma(division, args[[1]]*delta, args[[2]], args[[3]], args[[4]]*delta),
+## rstable=paste("rstable(division, ", args[2], ", ", args[3], ", ", args[4], ", ", args[5], ", ", args[6], ")")
+ rstable=rstable(division, args[[1]], args[[2]], args[[3]]*delta^(1/args[[1]]), args[[4]]*delta)
)
if(is.null(dZ)){ ##:: "otherwise"
cat(paste("Code \"", code, "\" not supported yet.\n", sep=""))
return(NULL)
}
- dZ <- eval(parse(text=dZ))
+ #print(parse(text=dZ))
+ #dZ <- eval(parse(text=dZ))
##:: calcurate difference eq.
for(i in 1:division){
Modified: pkg/yuima/R/yuima.model.R
===================================================================
--- pkg/yuima/R/yuima.model.R 2009-10-07 12:09:39 UTC (rev 4)
+++ pkg/yuima/R/yuima.model.R 2009-10-27 08:50:46 UTC (rev 5)
@@ -69,7 +69,7 @@
##::make type function list ####
CPlist <- c("dnorm", "dgamma", "dexp")
- codelist <- c("rIG", "rNIG", "rG", "rVG", "rS", "rTS", "rGH", "rgamma", "rbgamma", "rngamma", "rstable")
+ codelist <- c("rIG", "rNIG", "rgamma", "rbgamma", "rngamma", "rstable")
##::end make type function list ####
if(!length(measure.type)){
@@ -83,10 +83,10 @@
if( !length(jump.coeff) || !length(measure) ){
cat("\nmeasure type isn't matched with jump term.\n")
return(NULL)
- }else if(length(jump.coeff)!=1){
- cat("\nmulti dimentional jump term is not supported yet.\n")
- return(NULL)
- }
+ }#else if(length(jump.coeff)!=1){
+ # cat("\nmulti dimentional jump term is not supported yet.\n")
+ # return(NULL)
+ #}
if(measure.type=="CP"){ ##::CP
if(length(measure)!=2){
@@ -262,9 +262,14 @@
n.eqn2 <- n.eqn1
n.noise <- 1
}
+
+ ## TBC
+ n.eqn3 <- n.eqn1
+
if(!length(measure)){
n.eqn3 <- n.eqn1
}
+
if(n.eqn1 != n.eqn2 || n.eqn1 != n.eqn3){
cat("\nMalformed model, number of equations in the drift and diffusion do not match\n")
return(NULL)
Modified: pkg/yuima/man/adaBayes.Rd
===================================================================
--- pkg/yuima/man/adaBayes.Rd 2009-10-07 12:09:39 UTC (rev 4)
+++ pkg/yuima/man/adaBayes.Rd 2009-10-27 08:50:46 UTC (rev 5)
@@ -23,7 +23,7 @@
}
\author{YUIMA Project Team}
\note{
-This routine is not stable and accurate. Dr.Kamatani is now working for improvements.
+ This routine is not stable and accurate. Dr.Kamatani is now working for improvements.
}
\examples{
set.seed(1)
@@ -37,6 +37,9 @@
xinit=1,
true.parameter=c(0.6,0.3)
)
+
+ ## TBC: sample the path later
+
domain <- c(0,1)
## prior distribution of theta1
prior1 <- function(theta1,domain){
@@ -49,7 +52,7 @@
return( tmp )
}
- rand <- matrix(runif(10,0,1),10,2)
+ rand <- matrix(runif(20,0,1),10,2) ## initial points for MLE
QL <- 0
param <- numeric(2)
for( i in 1:5){
Modified: pkg/yuima/man/quasi-likelihood.Rd
===================================================================
--- pkg/yuima/man/quasi-likelihood.Rd 2009-10-07 12:09:39 UTC (rev 4)
+++ pkg/yuima/man/quasi-likelihood.Rd 2009-10-27 08:50:46 UTC (rev 5)
@@ -54,6 +54,40 @@
print("theta2 = .3, theta1 = .1")
print("ML estimator")
opt$par
+
+
+##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"))
+division <- 100
+ysamp <- setSampling(Terminal=(division)^(1/3), division=division)
+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=c(0.3, 0.1, 0.2, 0.1))
+
+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/((division)^(2/3)))
+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) )
+
+opt <- ml.ql(yuima, theta2, theta1, h=1/((division)^(2/3)), theta2.lim, theta1.lim)
+matrix(opt$par,2,2)
+
}
% Add one or more standard keywords, see file 'KEYWORDS' in the
Modified: pkg/yuima/man/simulate.Rd
===================================================================
--- pkg/yuima/man/simulate.Rd 2009-10-07 12:09:39 UTC (rev 4)
+++ pkg/yuima/man/simulate.Rd 2009-10-27 08:50:46 UTC (rev 5)
@@ -3,15 +3,16 @@
\title{Simulator function for multi-dimensional stochastic processes}
\description{Simulate multi-dimensional stochastic processes.}
\usage{
-simulate(yuima, xinit, true.parameter, space.discretized, increment)
+simulate(yuima, xinit, true.parameter, space.discretized, increment.W, increment.L)
}
\arguments{
\item{yuima}{an \code{yuima} object.}
\item{xinit}{initial value vector of state variables.}
- \item{true.parameter}{initial value vector of parameters.}
+ \item{true.parameter}{LIST of initial values of parameters.}
\item{space.discretized}{flag to switch to space-discretized Euler
Maruyama method.}
- \item{increment}{to specify increment for each time tics in advance.}
+ \item{increment.W}{to specify Wiener increment for each time tics in advance.}
+ \item{increment.L}{to specify Levy increment for each time tics in advance.}
}
\details{
\code{simulate} is a function to solve SDE using the Euler-Maruyama
@@ -42,6 +43,7 @@
# Solve SDEs using Euler-Maruyama method.
ou <- simulate(yuima=ou, xinit=1)
+x11()
plot(ou)
# A multi-dimensional (correlated) diffusion process.
@@ -119,9 +121,8 @@
yuima.mod <- simulate(yuima=yuima.mod,
xinit=1,
space.discretized=FALSE,
- increment=my.dW)
+ increment.W=my.dW)
if( !is.null(yuima.mod) ){
- x11()
plot(yuima.mod)
}
@@ -163,9 +164,8 @@
my.dW <- t(matrix(my.dW, nrow=division, ncol=yuima.obj at model@noise.number))
## Solve SDEs using Euler-Maruyama method.
-yuima.obj.path <- simulate(yuima=yuima.obj, space.discretized=FALSE, increment=my.dW)
+yuima.obj.path <- simulate(yuima=yuima.obj, space.discretized=FALSE, increment.W=my.dW)
if( !is.null(yuima.obj.path) ){
- x11()
plot(yuima.obj.path)
}
@@ -188,18 +188,30 @@
obj.sampling <- setSampling(Terminal=T, division=division)
obj.yuima <- setYuima(model=obj.model, sampling=obj.sampling)
-X <- simulate(yuima=obj.yuima, xinit=xinit, true.parameter=c(theta, sigma))
+X <- simulate(yuima=obj.yuima, xinit=xinit, true.parameter=list(theta, sigma))
plot(X)
##:: sample for Levy process ("code" type)
-## X_{t}= (sqrt(x+1) ) dt+ dw_t+ dZ_t
-obj.model <- setModel(drift="sqrt(x+1)", jump.coeff="1", measure.type="code", measure=list(df="rIG(z, 1, 0.1)"))
-obj.sampling <- setSampling(Terminal=1, division=5000)
+## dX_{t} = -x dt + dZ_t
+obj.model <- setModel(drift="-x", xinit=1, jump.coeff="1", measure.type="code", measure=list(df="rIG(z, 1, 0.1)"))
+obj.sampling <- setSampling(Terminal=10, division=10000)
obj.yuima <- setYuima(model=obj.model, sampling=obj.sampling)
result <- simulate(yuima=obj.yuima)
plot(result)
+##:: sample for multidimensional Levy process ("code" type)
+## dX = (theta - A X)dt + dZ,
+## theta=(theta_1, theta_2) = c(1,.5)
+## A=[a_ij], a_11 = 2, a_12 = 1, a_21 = 1, a_22=2
+ xinit <- c(1,1)
+ beta <- c(.1,.1)
+ Lambda <- matrix(c(1,0,0,1),2,2)
+ obj.model <- setModel(drift=c("1 - 2*x1-x2",".5-x1-2*x2"), xinit=c(1,1), solve.variable=c("x1","x2"), jump.coeff=c("1","1"), measure.type="code", measure=list(df="rNIG(z, alpha=1, beta=beta, delta=1,mu=c(0,0), Lambda=Lambda)"))
+ obj.sampling <- setSampling(Terminal=10, division=10000)
+ obj.yuima <- setYuima(model=obj.model, sampling=obj.sampling)
+ result <- simulate(yuima=obj.yuima,true.param=list(0,0,beta,Lambda))
+ plot(result)
}
More information about the Yuima-commits
mailing list