[Yuima-commits] r45 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Dec 16 11:08:27 CET 2009
Author: hinohide
Date: 2009-12-16 11:08:26 +0100 (Wed, 16 Dec 2009)
New Revision: 45
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/R/quasi-likelihood.R
pkg/yuima/man/quasi-likelihood.Rd
Log:
Extended the way to specify parameters for ml.ql function
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2009-12-12 14:14:14 UTC (rev 44)
+++ pkg/yuima/DESCRIPTION 2009-12-16 10:08:26 UTC (rev 45)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package
-Version: 0.0.76
-Date: 2009-12-12
+Version: 0.0.77
+Date: 2009-12-16
Depends: methods, zoo
Author: YUIMA Project Team.
Maintainer: Stefano M. Iacus <stefano.iacus at R-project.org>
Modified: pkg/yuima/R/quasi-likelihood.R
===================================================================
--- pkg/yuima/R/quasi-likelihood.R 2009-12-12 14:14:14 UTC (rev 44)
+++ pkg/yuima/R/quasi-likelihood.R 2009-12-16 10:08:26 UTC (rev 45)
@@ -123,36 +123,75 @@
##::theta1 : parameter in diffusion.
##::h : time width.
setGeneric("ql",
- function(yuima, theta2, theta1, h, print=FALSE)
+ function(yuima, theta2, theta1, h, print=FALSE, param)
standardGeneric("ql")
)
setMethod("ql", "yuima",
- function(yuima, theta2, theta1, h, print=FALSE){
+ function(yuima, theta2, theta1, h, print=FALSE, param){
##QLG <- ql.grad(yuima, theta2, theta1, h, print=FALSE)
##print(QLG)
-
if( missing(yuima)){
cat("\nyuima object is missing.\n")
return(NULL)
}
- if( missing(theta2) || missing(theta1)){
- cat("\nparameters of yuima.model is missing.\n")
- return(NULL)
+
+ ## param handling
+ if( missing(param) ){
+ if( missing(theta2) || missing(theta1) ){
+ cat("\nparameters of yuima.model is missing.\n")
+ return(NULL)
+ }
+ }else{
+ if( missing(theta2) && missing(theta1) ){
+ if( !is.list(param) ){
+ cat("\nparam must be list.\n")
+ return(NULL)
+ }
+
+ if( length(param)!=2){
+ cat("\nlength of param is strange.\n")
+ return(NULL)
+ }
+
+ ## get theta2 and theta1 from param
+ if( is.null(names(param)) ){
+ theta2 <- as.vector(param[[1]])
+ theta1 <- as.vector(param[[2]])
+ }
+ else if( sum( names(param)==c("theta2", "theta1") ) == 2 ){
+ theta2 <- as.vector(param[[1]])
+ theta1 <- as.vector(param[[2]])
+ }
+ else if( sum( names(param)==c("theta1", "theta2") ) == 2 ){
+ theta2 <- as.vector(param[[2]])
+ theta1 <- as.vector(param[[1]])
+ }
+ else{
+ cat("\nnames of param are strange.\n")
+ return(NULL)
+ }
+ }else{
+ cat("\nConflict in parameter specification method.\n")
+ return(NULL)
+ }
}
+ ## END param handling
+
if( missing(h)){
cat("\nlength of each time is missing.\n")
return(NULL)
}
+
if(length(yuima at model@parameter at drift)!=length(theta2)){
cat("\nlength of drift parameter is strange.\n")
return(NULL)
}
+
if(length(yuima at model@parameter at diffusion)!=length(theta1)){
cat("\nlength of diffusion parameter is strange.\n")
return(NULL)
}
-
d.size <- yuima at model@equation.number
division <- length(yuima)[1]
X <- as.matrix(onezoo(yuima))
@@ -266,19 +305,100 @@
##::ptheta1 : parameter in diffusion in the prevous mcmc step.
##::h : time width.
setGeneric("rql",
- function(yuima, theta2, theta1, ptheta2, ptheta1, h, print=FALSE)
+ function(yuima, theta2, theta1, ptheta2, ptheta1, h, print=FALSE, param, prevparam)
standardGeneric("rql")
)
setMethod("rql", "yuima",
- function(yuima, theta2, theta1, ptheta2, ptheta1, h, print=FALSE){
+ function(yuima, theta2, theta1, ptheta2, ptheta1, h, print=FALSE, param, prevparam){
if(missing(yuima)){
cat("\nyuima object is missing.\n")
return(NULL)
}
- if(missing(theta2) || missing(theta1)){
- cat("\nparameters of yuima.model is missing.\n")
- return(NULL)
+
+ ## param handling
+ if( missing(param) ){
+ if( missing(theta2) || missing(theta1) ){
+ cat("\nparameters of yuima.model is missing.\n")
+ return(NULL)
+ }
+ }else{
+ if( missing(theta2) && missing(theta1) ){
+ if( !is.list(param) ){
+ cat("\nparam must be list.\n")
+ return(NULL)
+ }
+
+ if( length(param)!=2){
+ cat("\nlength of param is strange.\n")
+ return(NULL)
+ }
+
+ ## get theta2 and theta1 from param
+ if( is.null(names(param)) ){
+ theta2 <- as.vector(param[[1]])
+ theta1 <- as.vector(param[[2]])
+ }
+ else if( sum( names(param)==c("theta2", "theta1") ) == 2 ){
+ theta2 <- as.vector(param[[1]])
+ theta1 <- as.vector(param[[2]])
+ }
+ else if( sum( names(param)==c("theta1", "theta2") ) == 2 ){
+ theta2 <- as.vector(param[[2]])
+ theta1 <- as.vector(param[[1]])
+ }
+ else{
+ cat("\nnames of param are strange.\n")
+ return(NULL)
+ }
+ }else{
+ cat("\nConflict in parameter specification method.\n")
+ return(NULL)
+ }
}
+ ## END param handling
+
+ ## prevparam handling
+ if( missing(prevparam) ){
+ if( missing(ptheta2) || missing(ptheta1) ){
+ cat("\nparameters of yuima.model is missing.\n")
+ return(NULL)
+ }
+ }else{
+ if( missing(ptheta2) && missing(ptheta1) ){
+ if( !is.list(prevparam) ){
+ cat("\nparam must be list.\n")
+ return(NULL)
+ }
+ if( length(prevparam)!=2){
+ cat("\nlength of param is strange.\n")
+ return(NULL)
+ }
+
+ ## get theta2 and theta1 from param
+ if( is.null(names(prevparam)) ){
+ ptheta2 <- as.vector(prevparam[[1]])
+ ptheta1 <- as.vector(prevparam[[2]])
+ }
+ else if( sum( names(prevparam)==c("ptheta2", "ptheta1") ) == 2 ){
+ ptheta2 <- as.vector(prevparam[[1]])
+ ptheta1 <- as.vector(prevparam[[2]])
+ }
+ else if( sum( names(prevparam)==c("ptheta1", "ptheta2") ) == 2 ){
+ ptheta2 <- as.vector(prevparam[[2]])
+ ptheta1 <- as.vector(prevparam[[1]])
+ }
+ else{
+ cat("\nnames of prevparam are strange.\n")
+ return(NULL)
+ }
+
+ }else{
+ cat("\nConflict in parameter specification method.")
+ return(NULL)
+ }
+ }
+ ## END prevparam handling
+
if(missing(h)){
cat("\nlength of each time is missing.\n")
return(NULL)
@@ -336,19 +456,59 @@
##::example: 0 <= theta1 <= 1 theta1.lim = matrix(c(0,1),1,2)
##::if theta1, theta2 are matrix, theta.lim can be defined matrix like rbind(c(0,1),c(0,1))
setGeneric("ml.ql",
- function(yuima, theta2, theta1, h, theta2.lim=matrix(c(0, 1), 1, 2), theta1.lim=matrix(c(0, 1), 1, 2), print=FALSE, BFGS=FALSE)
+ function(yuima, theta2, theta1, h, theta2.lim=matrix(c(0, 1), 1, 2), theta1.lim=matrix(c(0, 1), 1, 2), print=FALSE, BFGS=FALSE, param, interval)
standardGeneric("ml.ql")
)
setMethod("ml.ql", "yuima",
- function(yuima, theta2, theta1, h, theta2.lim=matrix(c(0, 1), 1, 2), theta1.lim=matrix(c(0, 1), 1, 2), print=FALSE, BFGS=FALSE){
+ function(yuima, theta2, theta1, h, theta2.lim=matrix(c(0, 1), 1, 2), theta1.lim=matrix(c(0, 1), 1, 2), print=FALSE, BFGS=FALSE, param, interval){
if( missing(yuima)){
cat("\nyuima object is missing.\n")
return(NULL)
}
- if( missing(theta2) || missing(theta1)){
- cat("\nparameters of yuima.model is missing.\n")
- return(NULL)
+
+ ## param handling
+ if( missing(param) ){
+ if( missing(theta2) || missing(theta1) ){
+ cat("\nparameters of yuima.model is missing.\n")
+ return(NULL)
+ }
+ }else{
+ if( missing(theta2) && missing(theta1) ){
+ if( !is.list(param) ){
+ cat("\nparam must be list.\n")
+ return(NULL)
+ }
+
+ if( length(param)!=2){
+ cat("\nlength of param is strange.\n")
+ return(NULL)
+ }
+
+ ## get theta2 and theta1 from param
+ if( is.null(names(param)) ){
+ theta2 <- as.vector(param[[1]])
+ theta1 <- as.vector(param[[2]])
+ }
+ else if( sum( names(param)==c("theta2", "theta1") ) == 2 ){
+ theta2 <- as.vector(param[[1]])
+ theta1 <- as.vector(param[[2]])
+ }
+ else if( sum( names(param)==c("theta1", "theta2") ) == 2 ){
+ theta2 <- as.vector(param[[2]])
+ theta1 <- as.vector(param[[1]])
+ }
+ else{
+ cat("\nnames of param are strange.\n")
+ return(NULL)
+ }
+ }else{
+ cat("\nConflict in parameter specification method.\n")
+ return(NULL)
+ }
}
+ ## END param handling
+
+
if(length(yuima at model@parameter at drift)!=length(theta2)){
cat("\nlength of drift parameter is strange.\n")
return(NULL)
@@ -361,10 +521,43 @@
cat("\nlength of each time is missing.\n")
return(NULL)
}
+
+ ## interval handling
+ if( !missing(interval) ){
+ if( missing(theta2.lim) && missing(theta2.lim) ){
+ if( !is.list(interval) ){
+ cat("\ninterval must be list.\n")
+ return(NULL)
+ }
+
+ theta2.len <- length(yuima at model@parameter at drift)
+ theta1.len <- length(yuima at model@parameter at diffusion)
+
+ if( length(interval) != (theta2.len+theta1.len) ){
+ cat("\nlength of interval is strange.\n")
+ return(NULL)
+ }
+
+ ## get theta2.lim and theta1.lim from interval
+ theta2.lim <- NULL
+ theta1.lim <- NULL
+ for(i in 1:theta2.len){
+ theta2.lim <- rbind(theta2.lim, interval[[i]])
+ }
+ for(i in 1:theta1.len){
+ theta1.lim <- rbind(theta1.lim, interval[[theta2.len+i]])
+ }
+ }else{
+ cat("\nConflict in parameter specification method.\n")
+ return(NULL)
+ }
+ }
+ ## END interval handling
+
if(is.matrix(theta2.lim) && is.matrix(theta1.lim)){
if(ncol(theta2.lim)!=2 || ncol(theta1.lim)!=2){
cat("\ntheta.lim is not available.\n")
- return(NULL)
+ return(NULL)
}
}
if( length(theta2)!=1 && length(theta2)!=nrow(theta2.lim)){
Modified: pkg/yuima/man/quasi-likelihood.Rd
===================================================================
--- pkg/yuima/man/quasi-likelihood.Rd 2009-12-12 14:14:14 UTC (rev 44)
+++ pkg/yuima/man/quasi-likelihood.Rd 2009-12-16 10:08:26 UTC (rev 45)
@@ -11,9 +11,9 @@
\description{Calculate the quasi-likelihood and estimate of the parameters of the
stochastic differential equation by the maximum likelihood method.}
\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,BFGS=FALSE)
-ql(yuima,theta2,theta1,h,print=FALSE)
-rql(yuima,theta2,theta1,ptheta2,ptheta1,h,print=FALSE)
+ml.ql(yuima,theta2,theta1,h,theta2.lim=matrix(c(0,1),1,2),theta1.lim=matrix(c(0,1),1,2),print=FALSE,BFGS=FALSE,param,interval)
+ql(yuima,theta2,theta1,h,print=FALSE,param)
+rql(yuima,theta2,theta1,ptheta2,ptheta1,h,print=FALSE,param,prevparam)
}
\arguments{
\item{yuima}{a yuima object.}
@@ -24,6 +24,9 @@
\item{ptheta2,ptheta1}{}
\item{print}{you can see a progress of the estimation when print is TRUE.}
\item{BFGS}{}
+ \item{param}{}
+ \item{interval}{}
+ \item{prevparam}{}
}
\details{
A function ql calculate the quasi-likelihood of a time series data X with any
@@ -51,6 +54,11 @@
QL <- ql(yuima, 0.8, 0.7, h=1/((division)^(2/3)))
QL
+## another way of parameter specification
+##param <- list(theta2=0.8, theta1=0.7)
+##QL <- ql(yuima, h=1/((division)^(2/3)), param=param)
+##QL
+
system.time(
opt <- ml.ql(yuima, 0.8, 0.7, h=1/((division)^(2/3)), c(0, 1), c(0, 1))
)
@@ -59,6 +67,16 @@
print("ML estimator")
opt$par
+## another way of parameter specification
+##interval <- list(theta2.lim=c(0,1), theta1.lim=c(0,1))
+##system.time(
+##opt <- ml.ql(yuima, h=1/((division)^(2/3)), param=param, interval=interval)
+##)
+##print("True param")
+##print("theta2 = .3, theta1 = .1")
+##print("ML estimator")
+##opt$par
+
system.time(
opt <- ml.ql(yuima, 0.8, 0.7, h=1/((division)^(2/3)), c(0, 1), c(0, 1), BFGS=TRUE)
)
@@ -67,7 +85,6 @@
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)
@@ -90,6 +107,11 @@
QL <- ql(yuima, theta2, theta1, h=1/((division)^(2/3)))
QL
+## another way of parameter specification
+#param <- list(theta2=theta2, theta1=theta1)
+#QL <- ql(yuima, h=1/((division)^(2/3)), param=param)
+#QL
+
theta2.1.lim <- c(0, 1)
theta2.2.lim <- c(0, 1)
theta1.1.lim <- c(0, 1)
@@ -102,6 +124,13 @@
)
opt$par
+## another way of parameter specification
+#interval <- list(theta2.1.lim, theta2.2.lim, theta1.1.lim, theta1.2.lim)
+#system.time(
+#opt <- ml.ql(yuima, h=1/((division)^(2/3)), param=param, interval=interval)
+#)
+#opt$par
+
system.time(
opt <- ml.ql(yuima, theta2, theta1, h=1/((division)^(2/3)), theta2.lim, theta1.lim, BFGS=TRUE)
)
More information about the Yuima-commits
mailing list