[Yuima-commits] r107 - in pkg/yuima: R man

[noreply at r-forge.r-project.org]

Author: iacus
Date: 2010-07-19 12:14:24 +0200 (Mon, 19 Jul 2010)
New Revision: 107

Added:
   pkg/yuima/R/lasso.R
   pkg/yuima/man/lasso.Rd
Log:
added lasso files as well

Added: pkg/yuima/R/lasso.R
===================================================================
--- pkg/yuima/R/lasso.R	                        (rev 0)
+++ pkg/yuima/R/lasso.R	2010-07-19 10:14:24 UTC (rev 107)
@@ -0,0 +1,61 @@
+# Initial version of lasso estimation for SDEs
+
+lasso <- function(yuima, lambda0, start, ...){
+	
+	call <- match.call()
+	
+	if( missing(yuima))
+		yuima.stop("yuima object 'yuima' is missing.")
+	
+	if( missing(lambda0) ){
+	 pars <- yuima at model@parameter at all	
+	 lambda0 <- rep(1, length(pars))
+	 names(lambda0) <- pars
+	 lambda0 <- as.list(lambda0)	
+	}
+	
+## 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.")
+	
+	fail <- lapply(lambda0, function(x) as.numeric(NA))
+	
+	cat("\nLooking for MLE estimates...\n")
+	fit <- try(qmle(yuima, start=start,...), silent=TRUE)
+	if(class(fit)=="try-error")
+	return(list(mle=fail, sd.mle=NA, lasso=fail, sd.lasso=NA))
+	
+	SIGMA <- try( sqrt(diag(vcov(fit))), silent=TRUE)
+	if(class(SIGMA)=="try-error")
+	return(list(mle=fail, sd.mle=NA, lasso=fail, sd.lasso=NA))
+	
+	
+	theta.mle <- coef(fit)
+	
+	H <- try( solve(vcov(fit)), silent=TRUE) 
+	
+	if(class(H)=="try-error")
+	return(list(mle=fail, sd.mle=NA, lasso=fail, sd.lasso=NA))
+	
+	
+	lambda <- unlist(lambda0[names(theta.mle)])/abs(theta.mle)
+	
+	f2 <- function( theta ) as.numeric( t(theta-theta.mle) %*% H %*% (theta-theta.mle) + lambda %*% abs(theta) )
+	
+	cat("\nPerforming LASSO estimation...\n")
+	
+	fit2 <- try( optim(theta.mle, f2, hessian=TRUE,..., 
+					   control = list(maxit=30000, temp=2000, REPORT=500)), silent=TRUE) 
+	if(class(fit2)=="try-error")
+	return(list(mle=fail, sd.mle=NA, lasso=fail, sd.lasso=NA))
+	
+	theta.lasso <- fit2$par
+	
+	SIGMA1 <- try(sqrt(diag(solve(fit2$hessian))), silent=TRUE)
+	
+	if(class(SIGMA1)=="try-error")
+	return(list(mle=fail, sd.mle=NA, lasso=fail, sd.lasso=NA))
+	
+	return(list(mle=theta.mle, sd.mle=SIGMA, lasso=theta.lasso, sd.lasso=SIGMA1,call=call, lambda0=lambda0))
+}

Added: pkg/yuima/man/lasso.Rd
===================================================================
--- pkg/yuima/man/lasso.Rd	                        (rev 0)
+++ pkg/yuima/man/lasso.Rd	2010-07-19 10:14:24 UTC (rev 107)
@@ -0,0 +1,63 @@
+\name{lasso}
+\alias{lasso}
+\title{Adaptive LASSO estimation for stochastic differential equations}
+\description{Adaptive LASSO estimation for stochastic differential equations.}
+\usage{
+lasso(yuima, lambda0, start, ...)
+}
+\arguments{
+  \item{yuima}{a yuima object.}
+  \item{lambda0}{a named list with penalty for each parameter.}
+  \item{start}{initial values to be passed to the optimizer.}
+  \item{...}{passed to \code{\link{optim}} method. See Examples.}
+}
+\details{
+  
+  \code{lasso} behaves more likely the standard \code{\link{qmle}} function in  and
+  argument \code{method} is one of the methods available in \code{\link{optim}}.
+
+  From initial guess of QML estimates, performs adaptive LASSO estimation using 
+  the Least Squares Approximation (LSA) as in Wang and Leng (2007, JASA).
+  
+  
+}
+\value{
+  \item{ans}{a list with both QMLE and LASSO estimates.}
+}
+\author{The YUIMA Project Team}
+\examples{
+##multidimension case
+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)
+
+truep <- list(theta1.1=0.6, theta1.2=0,theta2.1=0.5, theta2.2=0)
+yuima <- simulate(yuima, xinit=c(1, 1), 
+ true.parameter=truep)
+
+
+est <- lasso(yuima, start=list(theta2.1=0.8, theta2.2=0.2, theta1.1=0.7, theta1.2=0.1),
+ lower=list(theta1.1=1e-10,theta1.2=1e-10,theta2.1=.1,theta2.2=1e-10),
+ upper=list(theta1.1=4,theta1.2=4,theta2.1=4,theta2.2=4), method="L-BFGS-B")
+
+# TRUE
+unlist(truep)
+
+# QMLE
+round(est$mle,3)
+
+# LASSO
+round(est$lasso,3) 
+}
+ 
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
+\keyword{ts}