[Yuima-commits] r602 - in pkg/yuima: . R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Apr 12 13:19:30 CEST 2017
Author: kyuta
Date: 2017-04-12 13:19:28 +0200 (Wed, 12 Apr 2017)
New Revision: 602
Added:
pkg/yuima/src/euler.c
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NEWS
pkg/yuima/R/sim.euler.R
Log:
fix a bug in sim.euler.R
(re-)added euler.c
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2017-03-27 05:56:42 UTC (rev 601)
+++ pkg/yuima/DESCRIPTION 2017-04-12 11:19:28 UTC (rev 602)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project Package for SDEs
-Version: 1.6.1
+Version: 1.6.2
Depends: R(>= 2.10.0), methods, zoo, stats4, utils, expm, cubature, mvtnorm
Imports: Rcpp (>= 0.12.1)
Author: YUIMA Project Team
Modified: pkg/yuima/NEWS
===================================================================
--- pkg/yuima/NEWS 2017-03-27 05:56:42 UTC (rev 601)
+++ pkg/yuima/NEWS 2017-04-12 11:19:28 UTC (rev 602)
@@ -45,5 +45,7 @@
2016/12/16: added rGIG, rGH, dGIG and dGH in rng.R and the corresponding c language file YU
2017/01/25: modified sim.euler.R and added euler.c to implement the Euler-Maruyama scheme by the C code (only the diffusion case)
2017/02/23: modified sim.euler.R and removed euler.c due to a memory corruption
- fix a bug in sim.euler.R
-2017/03/27: added IC.R and qmleLevy.R
\ No newline at end of file
+ fix a bug in sim.euler.R
+2017/03/27: added IC.R and qmleLevy.R
+2017/04/12: fix a bug in sim.euler.R
+ (re-)added euler.c
\ No newline at end of file
Modified: pkg/yuima/R/sim.euler.R
===================================================================
--- pkg/yuima/R/sim.euler.R 2017-03-27 05:56:42 UTC (rev 601)
+++ pkg/yuima/R/sim.euler.R 2017-04-12 11:19:28 UTC (rev 602)
@@ -55,7 +55,9 @@
##:: set time step
# delta <- Terminal/n
delta <- yuima at sampling@delta
-
+
+ ## moving the following lines enclosed if(0) after Levy (YK, Apr 12, 2017)
+ if(0){
##:: 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)))
@@ -93,16 +95,16 @@
X_mat <- matrix(0, d.size, n+1)
X_mat[,1] <- dX
-
+
if(has.drift){ ##:: consider drift term to be one of the diffusion term(dW=1)
dW <- rbind( rep(1, n)*delta , dW)
}
+ }## if(0) finish
-
if(!length(sdeModel at measure.type)){ ##:: Wiener Proc
##:: using Euler-Maruyama method
- #if(0){ # old version (Jan 25, 2017)
+ if(0){ # old version (Jan 25, 2017)
if(var.in.diff & (!is.Poisson(sdeModel))){ ##:: diffusions have state variables and it is not Poisson
##:: calcurate difference eq.
for( i in 1:n){
@@ -155,9 +157,9 @@
X_mat[, i+1] <- dX
}
}
- #}
+ }
- if(0){ # currently ignored due to a bug (YK, Feb 23, 2017)
+ #if(0){ # currently ignored due to a bug (YK, Feb 23, 2017)
# new version (Jan 25, 2017)
b <- parse(text=paste("c(",paste(as.character(V0),collapse=","),")"))
vecV <- parse(text=paste("c(",paste(as.character(unlist(V)),collapse=","),")"))
@@ -165,10 +167,53 @@
X_mat <- .Call("euler", dX, Initial, as.integer(r.size),
rep(1, n) * delta, dW, modeltime, modelstate, quote(eval(b, env)),
quote(eval(vecV, env)), env, new.env())
- }
+ #}
#tsX <- ts(data=t(X_mat), deltat=delta , start=0)
tsX <- ts(data=t(X_mat), deltat=delta , start = yuima at sampling@Initial) #LM
}else{ ##:: Levy
+ ### add (YK, Apr 12, 2017)
+ ##:: 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)))
+ #print(is.Poisson(sdeModel))
+
+ ##:: function to calculate coefficients of dW(including drift term)
+ ##:: common used in Wiener and CP
+ p.b <- function(t, X=numeric(d.size)){
+ ##:: assign names of variables
+ for(i in 1:length(modelstate)){
+ assign(modelstate[i], X[i], env)
+ }
+ assign(modeltime, t, env)
+ ##:: solve diffusion term
+ if(has.drift){
+ tmp <- matrix(0, d.size, r.size+1)
+ for(i in 1:d.size){
+ tmp[i,1] <- eval(V0[i], env)
+ for(j in 1:r.size){
+ tmp[i,j+1] <- eval(V[[i]][j],env)
+ }
+ }
+ } else { ##:: no drift term (faster)
+ tmp <- matrix(0, d.size, r.size)
+ if(!is.Poisson(sdeModel)){ # we do not need to evaluate diffusion
+ for(i in 1:d.size){
+ for(j in 1:r.size){
+ tmp[i,j] <- eval(V[[i]][j],env)
+ } # for j
+ } # foh i
+ } # !is.Poisson
+ } # else
+ return(tmp)
+ }
+
+ X_mat <- matrix(0, d.size, n+1)
+ X_mat[,1] <- dX
+
+ if(has.drift){ ##:: consider drift term to be one of the diffusion term(dW=1)
+ dW <- rbind( rep(1, n)*delta , dW)
+ }
+
JP <- sdeModel at jump.coeff
mu.size <- length(JP)
# cat("\n Levy\n")
Added: pkg/yuima/src/euler.c
===================================================================
--- pkg/yuima/src/euler.c (rev 0)
+++ pkg/yuima/src/euler.c 2017-04-12 11:19:28 UTC (rev 602)
@@ -0,0 +1,89 @@
+//
+// euler.c
+//
+//
+// Created by Yuta Koike on 2016/01/23.
+//
+//
+
+#include <R.h>
+#include <Rmath.h>
+#include <Rdefines.h>
+#include <Rinternals.h>
+
+SEXP euler(SEXP x0, SEXP t0, SEXP R, SEXP dt, SEXP dW, SEXP modeltime, SEXP modelstate,
+ SEXP drift, SEXP diffusion, SEXP env, SEXP rho){
+
+ int i, j, k, n, d, r;
+ double *rdt, *rdW, *rX, *rx0, *b, *sigma;
+ SEXP X, xpar, tpar;
+ //SEXP X, xpar, tpar, xvar, tvar;
+
+ PROTECT(x0 = AS_NUMERIC(x0));
+ rx0 = REAL(x0);
+ PROTECT(dt = AS_NUMERIC(dt));
+ rdt = REAL(dt);
+ PROTECT(dW = AS_NUMERIC(dW));
+ rdW = REAL(dW);
+
+ d = length(x0);
+ n = length(dt);
+
+ r = *INTEGER(R);
+
+ /* PROTECT(X = allocVector(REALSXP, n + 1));
+ rX = REAL(X);
+ rX[0] = *REAL(x0);*/
+ PROTECT(X = allocMatrix(REALSXP, d, n + 1));
+ rX = REAL(X);
+
+ for (j = 0; j < d; j++) {
+ rX[j] = rx0[j];
+ }
+
+ PROTECT(tpar = allocVector(REALSXP, 1));
+ //PROTECT(xpar = allocVector(REALSXP, 1));
+
+ PROTECT(t0 = AS_NUMERIC(t0));
+ REAL(tpar)[0] = REAL(t0)[0]; /* initial time */
+
+ for (i = 0; i < n; i++) {
+
+ /* assign the current variables */
+ defineVar(installChar(STRING_ELT(modeltime, 0)), tpar, env);
+ // PROTECT(tvar = STRING_ELT(modeltime, 0));
+ //defineVar(installChar(tvar), tpar, env);
+
+ for (j = 0; j < d; j++) {
+ PROTECT(xpar = allocVector(REALSXP, 1));
+ REAL(xpar)[0] = rX[j + i * d];
+ //defineVar(installChar(STRING_ELT(modelstate, j)), duplicate(xpar), env);
+ //PROTECT(xvar = STRING_ELT(modelstate, j));
+ //defineVar(installChar(xvar), duplicate(xpar), env);
+ defineVar(installChar(STRING_ELT(modelstate, j)), xpar, env);
+ UNPROTECT(1);
+ }
+
+ /*defineVar(install("env"), env, rho);*/
+
+ /* evaluate coefficients */
+ b = REAL(eval(drift, rho));
+ sigma = REAL(eval(diffusion, rho));
+
+ for (j = 0; j < d; j++) {
+ rX[j + (i + 1) * d] = rX[j + i * d] + b[j] * rdt[i];
+ for (k = 0; k < r; k++) {
+ rX[j + (i + 1) * d] += sigma[k + j * r] * rdW[k + i * r];
+ }
+ }
+
+ /*defineVar(install("t"), tpar, rho);*/
+ /*rX[i + 1] = rX[i] + *REAL(eval(drift, rho)) * rdt[i] + *REAL(eval(diffusion, rho)) * rdW[i];*/
+ /*rX[i + 1] = rX[i] + *REAL(eval(drift, rho)) * REAL(dt)[i] + *REAL(eval(diffusion, rho)) * REAL(dW)[i];*/
+
+ REAL(tpar)[0] += rdt[i];
+ }
+
+ UNPROTECT(6);
+ return(X);
+}
More information about the Yuima-commits
mailing list