[Yuima-commits] r343 - in pkg/yuima: . R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Oct 23 13:37:46 CEST 2014
Author: lorenzo
Date: 2014-10-23 13:37:45 +0200 (Thu, 23 Oct 2014)
New Revision: 343
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/R/qmle.R
pkg/yuima/src/carmafilter.c
Log:
update qmle.R for carma model
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2014-10-20 08:21:58 UTC (rev 342)
+++ pkg/yuima/DESCRIPTION 2014-10-23 11:37:45 UTC (rev 343)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package for SDEs
-Version: 1.0.41
-Date: 2014-10-18
+Version: 1.0.42
+Date: 2014-10-23
Depends: methods, zoo, stats4, utils, expm, cubature, mvtnorm
Author: YUIMA Project Team
Maintainer: Stefano M. Iacus <stefano.iacus at unimi.it>
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2014-10-20 08:21:58 UTC (rev 342)
+++ pkg/yuima/R/qmle.R 2014-10-23 11:37:45 UTC (rev 343)
@@ -1823,14 +1823,10 @@
#carma.kalman<-function(y, tt, p, q, a,bvector, sigma){
-carma.kalman<-function(y, u, p, q, a,bvector, sigma, times.obs, V_inf0){
-
- # V_inf0<-matrix(diag(rep(1,p)),p,p)
-
+carma.kalman<-function(y, u, p, q, a,bvector, sigma, times.obs, V_inf0){
+ # V_inf0<-matrix(diag(rep(1,p)),p,p)
A<-MatrixA(a)
# u<-diff(tt)[1]
-
-
# Amatx<-yuima.ca2rma.eigen(A)
# expA<-Amatx$vectors%*%expm(diag(Amatx$values*u),
# method="Pade",
@@ -1844,59 +1840,66 @@
expA<-expm(A*u,method="Pade",order=6, trySym=FALSE, do.sparseMsg = FALSE)
# }
#expA<-yuima.exp(A*u)
-
- v<-as.numeric(V_inf0[upper.tri(V_inf0,diag=TRUE)])
-
- ATrans<-t(A)
- matrixV<-matrix(0,p,p)
- #l.dummy<-c(rep(0,p-1),1)
- l<-rbind(matrix(rep(0,p-1),p-1,1),1)
- #l<-matrix(l.dummy,p,1)
- #lTrans<-matrix(l.dummy,1,p)
- lTrans<-t(l)
- elForVInf<-new.env()
- elForVInf$A<-A
- elForVInf$ATrans<-ATrans
- elForVInf$lTrans<-lTrans
- elForVInf$l<-l
- elForVInf$matrixV<-matrixV
-
- elForVInf$sigma<-sigma
-# elForVInf<-list(A=A,
-# ATrans=ATrans,
-# lTrans=lTrans,
-# l=l,
-# matrixV=matrixV,
-# sigma=sigma)
+# v<-as.numeric(V_inf0[upper.tri(V_inf0,diag=TRUE)])
+#
+# ATrans<-t(A)
+# matrixV<-matrix(0,p,p)
+# #l.dummy<-c(rep(0,p-1),1)
+# l<-rbind(matrix(rep(0,p-1),p-1,1),1)
+# #l<-matrix(l.dummy,p,1)
+# #lTrans<-matrix(l.dummy,1,p)
+# lTrans<-t(l)
+# elForVInf<-new.env()
+# elForVInf$A<-A
+# elForVInf$ATrans<-ATrans
+# elForVInf$lTrans<-lTrans
+# elForVInf$l<-l
+# elForVInf$matrixV<-matrixV
+#
+# elForVInf$sigma<-sigma
+# # elForVInf<-list(A=A,
+# # ATrans=ATrans,
+# # lTrans=lTrans,
+# # l=l,
+# # matrixV=matrixV,
+# # sigma=sigma)
+# #
+#
+# # V_inf_vect<-nlm(yuima.Vinfinity, v,
+# # elForVInf = elForVInf)$estimate
+#
+# V_inf_vect<-nlminb(start=v,objective=yuima.Vinfinity, elForVInf = elForVInf)$par
+ # V_inf_vect<-optim(par=v,fn=yuima.Vinfinity,method="L-BFGS-B", elForVInf = elForVInf)$par
+# V_inf<-diag(p)
+#
+# V_inf[upper.tri(V_inf,diag=TRUE)]<-V_inf_vect
+# V_inf<-as.matrix(forceSymmetric(V_inf))
#
- V_inf_vect<-nlm(yuima.Vinfinity, v,
- elForVInf = elForVInf)$estimate
-# V_inf_vect<-nlminb(start=v,objective=yuima.Vinfinity, elForVInf = elForVInf)$par
-# V_inf_vect<-optim(par=v,fn=yuima.Vinfinity,method="L-BFGS-B", elForVInf = elForVInf)$par
- V_inf<-matrix(0,p,p)
+# V_inf[abs(V_inf)<=1.e-06]=0
+# pp = as.integer(length(a))
+# vinf = matrix(ncol = pp, nrow = pp, 0)
+# xx = .Fortran("vinfinity", a, pp, vinf, PACKAGE = "yuima")
+# We need to write this part using Yuima
+# V_inf<-sigma^2 * xx[[3]]
+
+ V_inf<-V0inf(a,p)*sigma^2
- V_inf[upper.tri(V_inf,diag=TRUE)]<-V_inf_vect
- V_inf<-as.matrix(forceSymmetric(V_inf))
-# V_inf[lower.tri(V_inf)]<-V_inf[upper.tri(V_inf)]
-
- V_inf[abs(V_inf)<=1.e-06]=0
-
expAT<-t(expA)
#SIGMA_err<-V_inf-expA%*%V_inf%*%t(expA)
# input: V_inf, p, expA
# output: SigMatr (pre-alloc in R)
#
- SigMatr <- .Call("carma_tmp", V_inf, as.integer(p), expA, PACKAGE="yuima")
+# SigMatr <- .Call("carma_tmp", V_inf, as.integer(p), expA, PACKAGE="yuima")
#
-# SigMatr <- V_inf - expA %*% V_inf %*% expAT
+ SigMatr <- V_inf - expA %*% V_inf %*% expAT
# cat("\nSigMatr\n")
# print(SigMatr)
# stop("")
#
- statevar<-matrix(rep(0, p),p,1)
+ statevar<-numeric(length=p)
Qmatr<-SigMatr+0
# set
@@ -1909,20 +1912,25 @@
#SigMatr <- V_inf
- zc <- matrix(bvector,1,p)
- loglstar <- 0
+ zc <- matrix(bvector,1,p)
+ loglstar <- 0
-# zcT<-matrix(bvector,p,1)
- zcT <- t(zc)
+ zcT<-matrix(bvector,p,1)
+ #zcT <- t(zc)
# Cycle_Carma(SEXP StateVar, SEXP ExpA, SEXP Times.Obs, SEXP P,
# SEXP Qmatr, SEXP SigMatr, SEXP Zc, SEXP Logstar)
- # sd_2<-0
-
- loglstar<- .Call("Cycle_Carma", y, statevar, expA, as.integer(times.obs),
- as.integer(p), Qmatr, SigMatr, zc,
- PACKAGE="yuima")
- return(list(loglstar=as.numeric(loglstar[1,1])-0.5*log(2*pi)*times.obs,s2hat=loglstar[2,1]))
+ # sd_2<-0
+ Result<-numeric(length=2)
+ Kgain<-numeric(length=p)
+ dum_zc<-numeric(length=p)
+ Mat22int<-numeric(length=(p*p))
+# #
+ loglstar<- .Call("Cycle_Carma", y, statevar, expA, as.integer(length(y)),
+ as.integer(p), Qmatr, SigMatr, bvector, Result, Kgain,
+ dum_zc, Mat22int,
+ PACKAGE="yuima")
+ return(list(loglstar=loglstar[1]-0.5*log(2*pi)*times.obs,s2hat=loglstar[2]))
# statevar, expA, times.obs, Qmatr, SigMatr, zc, logstar
@@ -1933,20 +1941,53 @@
# statevar <- expA %*% statevar
# SigMatr <- expA %*% SigMatr %*% expAT + Qmatr
# # forecast
-# Uobs <- y[t] - zc %*% statevar
+# #Uobs <- y[t] - zc %*% statevar
+# Uobs <- y[t] - sum(bvector * statevar)
# dum.zc <- zc %*% SigMatr
-# sd_2 <- dum.zc %*% zcT
+# sd_2 <- sum(dum.zc * bvector)
# Inv_sd_2 <- 1/sd_2
-# #correction
-# Kgain <- SigMatr %*% zcT %*% Inv_sd_2
-# statevar <- statevar+Kgain %*% Uobs
+# #correction Kgain <- SigMatr %*% zcT %*% Inv_sd_2
+# Kgain <- SigMatr %*% zcT *Inv_sd_2
+# statevar <- statevar+Kgain * Uobs
# SigMatr <- SigMatr - Kgain %*% dum.zc
-# term_int<- -0.5 * (log(sd_2)+ Uobs %*% Uobs %*% Inv_sd_2)
+# term_int<- -0.5 * (log(sd_2)+ Uobs * Uobs *Inv_sd_2)
# loglstar <- loglstar + term_int
# }
# return(list(loglstar=loglstar-0.5*log(2*pi)*times.obs,s2hat=sd_2))
}
+V0inf<-function(a,p){
+ # This code is based on the paper A continuous-time ARMA process Tsai-Chan 2000
+ # we need to find the values along the diagonal
+ #l<-c(numeric(length=(p-1)),0.5)
+ # B_{p*p}V^{*}_{p*1}=-sigma^2*l/2
+ B<-matrix(0,nrow=p,ncol=p)
+ aa <- -rev(a)
+ for(i in 1:p){
+ # Condition on B
+ for(j in 1:p){
+ if ((2*j-i) %in% c(1:p)){
+ B[i,j]<-(-1)^(j-i)*aa[2*j-i]
+ }
+ if((2*j-i)==(p+1)){
+ B[i,j]<-(-1)^(j-i-1)
+ }
+ }
+ }
+ Vdiag <- -solve(B)[,p]*0.5
+ V <- diag(Vdiag)
+ # we insert the values outside the diagonal
+ for(i in 1:p){
+ for(j in (i+1):p){
+ if((i+j) %% 2 == 0){ # if even
+ V[i,j]=(-1)^((i-j)/2)*V[(i+j)/2,(i+j)/2]
+ V[j,i]=V[i,j]
+ }
+ }
+ }
+ return(V)
+}
+
# CycleCarma<-function(y, statevar, expA, times.obs=integer(),
# p=integer(), Qmatr, SigMatr, zc, loglstar){
# # expAT=t(expA)
Modified: pkg/yuima/src/carmafilter.c
===================================================================
--- pkg/yuima/src/carmafilter.c 2014-10-20 08:21:58 UTC (rev 342)
+++ pkg/yuima/src/carmafilter.c 2014-10-23 11:37:45 UTC (rev 343)
@@ -44,7 +44,8 @@
SEXP carma_tmp(SEXP V, SEXP P, SEXP A);
SEXP Cycle_Carma(SEXP Y, SEXP StateVar, SEXP ExpA, SEXP Times_Obs, SEXP P,
- SEXP Qmatr, SEXP SigMatr, SEXP Zc);
+ SEXP Qmatr, SEXP SigMatr, SEXP Zc, SEXP Result,SEXP Kgain,
+ SEXP dum_zc, SEXP Mat22int);
@@ -105,37 +106,36 @@
}
SEXP Cycle_Carma(SEXP Y, SEXP StateVar, SEXP ExpA, SEXP Times_Obs, SEXP P,
- SEXP Qmatr, SEXP SigMatr, SEXP Zc){
+ SEXP Qmatr, SEXP SigMatr, SEXP Zc, SEXP Result, SEXP Kgain,
+ SEXP dum_zc, SEXP Mat22int){
/* Declare Integer Variable */
int times_obs, p;
- int i, j, h, k, t;
+ int i, j, h, t;
+
/* Declare pointer */
double *rY, *rStateVar, *rExpA, *rQmatr, *rSigMatr, *rZc;
- double *rKgain, *rdum_zc, *rMat21, *rMat22int, *rMat22est, *rResult;
- double ppp=0;
+ double *rKgain, *rdum_zc, *rMat22int, *rResult;
double Uobs=0;
- double sd_2=0;
- double LogL=0;
- double dummy=0;
-
+ double dummy=0;
+ double dummy1=0;
/* Declare SEXP */
- SEXP Kgain, dum_zc, Mat21, Mat22int, Mat22est, Result;
+
/* Check the type of variables*/
- if(!isInteger(P)) error("`P' must be integer");
- if(!isInteger(Times_Obs)) error("`Times_Obs' must be integer");
- if(!isNumeric(Y)) error("`Y' must be numeric");
- if(!isNumeric(StateVar)) error("`StateVar' must be numeric");
- if(!isNumeric(ExpA)) error("`ExpA' must be numeric");
- if(!isNumeric(Qmatr)) error("`Qmatr' must be numeric");
- if(!isNumeric(SigMatr)) error("`SigMatr' must be numeric");
- if(!isNumeric(Zc)) error("`Zc' must be numeric");
+ // if(!isInteger(P)) error("`P' must be integer");
+ // if(!isInteger(Times_Obs)) error("`Times_Obs' must be integer");
+ // if(!isNumeric(Y)) error("`Y' must be numeric");
+ // if(!isNumeric(StateVar)) error("`StateVar' must be numeric");
+ // if(!isNumeric(ExpA)) error("`ExpA' must be numeric");
+ // if(!isNumeric(Qmatr)) error("`Qmatr' must be numeric");
+ // if(!isNumeric(SigMatr)) error("`SigMatr' must be numeric");
+ // if(!isNumeric(Zc)) error("`Zc' must be numeric");
/* Protect Objects */
@@ -157,7 +157,20 @@
PROTECT(Zc = AS_NUMERIC(Zc));
rZc = REAL(Zc);
+
+ PROTECT(Result = AS_NUMERIC(Result));
+ rResult = REAL(Result);
+
+ PROTECT(Kgain = AS_NUMERIC(Kgain));
+ rKgain = REAL(Kgain);
+
+ PROTECT(dum_zc = AS_NUMERIC(dum_zc));
+ rdum_zc = REAL(dum_zc);
+
+ PROTECT(Mat22int = AS_NUMERIC(Mat22int));
+ rMat22int = REAL(Mat22int);
+
/*PROTECT(Logstar = AS_NUMERIC(Logstar));
rLoglstar = REAL(Logstar);*/
@@ -166,27 +179,9 @@
times_obs = *INTEGER(Times_Obs);
p = *INTEGER(P);
- // PROTECT(Mat21 = allocMatrix(REALSXP, p, 1));
- // rMat21 = REAL(Mat21);
- PROTECT(Kgain = allocMatrix(REALSXP, p, 1));
- rKgain = REAL(Kgain);
-
- PROTECT(dum_zc = allocMatrix(REALSXP, p,1));
- rdum_zc = REAL(dum_zc);
-
- PROTECT(Mat22int = allocMatrix(REALSXP, p, p));
- rMat22int = REAL(Mat22int);
-
- PROTECT(Mat22est = allocMatrix(REALSXP, p, p));
- rMat22est = REAL(Mat22est);
-
- PROTECT(Result = allocMatrix(REALSXP, 2, 1));
- rResult = REAL(Result);
-
-
- /* Main Code */
- /* Dimension of Inputs:
+ /* Main Code
+ Dimension of Inputs:
Y = Vector p dimension;
StateVar = matrix p x 1;
ExpA = matrix p x p;
@@ -200,9 +195,6 @@
for(t in 1:rtime_obs){*/
rResult[0]=0;
- /* for(i=0; i<p; i++){
- rdum_zc[i]=0;
- }*/
for(t=0; t<times_obs; t++){
/*t=0;*/
@@ -214,9 +206,7 @@
dummy += rExpA[i+j*p]*rStateVar[j];
}
rStateVar[i] = dummy;
- /* rMat21[i]=0; statevar <- expA %*% statevar */
- /*}
- SigMatr <- expA %*% SigMatr %*% expAT + Qmatr: pxp = pxp pxp pxp
+ /*SigMatr <- expA %*% SigMatr %*% expAT + Qmatr: pxp = pxp pxp pxp
First We compute rMat22int <- expA %*% SigMatr : pxp = pxp pxp
for(i=0; i<p; i++){*/
for(j=0; j<p; j++){
@@ -229,223 +219,59 @@
}
/* Second We compute rMat22est <- rMat22int %*% t(expA)
+ Qmatr: pxp = pxp pxp + pxp */
+ Uobs = 0;
for(i=0; i<p; i++){
for(j=0; j<p; j++){
- rMat22est[i+j*p]=0;
+ rSigMatr[i+j*p]=0;
for(h=0; h<p; h++){
- rMat22est[i+j*p]=rMat22est[i+j*p]+rMat22int[i+h*p]*rExpA[j+h*p];
+ rSigMatr[i+j*p] += rMat22int[i+h*p]*rExpA[j+h*p];
}
- rSigMatr[i+j*p]=rMat22est[i+j*p]+rQmatr[i+j*p];
+ rSigMatr[i+j*p] = rSigMatr[i+j*p]+rQmatr[i+j*p];
}
-
+ /*forecast
+ Uobs <- y[t] - zc %*% statevar # 1 - 1Xp px*/
+ Uobs += rZc[i]*rStateVar[i];
}
-
- /*forecast*/
-
- /*Uobs <- y[t] - zc %*% statevar # 1 - 1Xp px1*/
- rMat22est[0]=0;
- for(i=0; i<p; i++){
- rMat22est[0]=rMat22est[0]+rZc[i]*rStateVar[i];
- }
- Uobs=rY[t]-rMat22est[0];
- /*rMat22est[0]=0;*/
+ Uobs = rY[t]-Uobs;
/* dum.zc <- zc %*% SigMatr # 1xp pxp */
+ rResult[1] =0;
for(i=0; i<p; i++){
dummy = 0;
for(j=0; j<p; j++){
dummy += rSigMatr[i+j*p]*rZc[j];
}
rdum_zc[i]=dummy;
- }
-
-
- /*rResult[1]=rSigMatr[0];*/
/* Sd_2 <- dum.zc %*% zcT # 1xp px1 */
- rResult[1] =0;
- for(i=0; i<p; i++){
- rResult[1] = rResult[1]+rdum_zc[i]*rZc[i];
+ rResult[1] += rdum_zc[i]*rZc[i];
}
-
/* correction */
/* Kgain <- SigMatr %*% zcT %*% 1/Sd_2 # pxp px1*1 */
for(i=0; i<p; i++){
- dummy=0;
+ dummy1=0;
for(j=0; j<p; j++){
- dummy += rSigMatr[i+j*p]*rZc[j];
+ dummy1 += rSigMatr[i+j*p]*rZc[j];
}
- rKgain[i]=dummy/rResult[1];
- /*rMat21[i]=0;*/
- }
+ rKgain[i]=dummy1/rResult[1];
/* statevar <- statevar+Kgain %*% Uobs # px1+px1 */
- for(i=0; i<p; i++){
- rStateVar[i] = rStateVar[i] + rKgain[i]*Uobs;
- }
-
-
- /* SigMatr <- SigMatr - Kgain %*% dum.zc # pxp-px1 1xp */
- for(i=0; i<p; i++){
+ rStateVar[i] += rKgain[i]*Uobs;
+ /*SigMatr <- SigMatr - Kgain %*% dum.zc # pxp-px1 1xp*/
for(j=0; j<p; j++){
- rSigMatr[i+j*p]=rSigMatr[i+j*p]-rKgain[i]*rdum_zc[j];
+ rSigMatr[i+j*p] += -rKgain[i]*rdum_zc[j];
}
}
- /*for(i=0; i<p*p; i++){
- rMat22int[i];
- }*/
/*term_int = -0.5 * (log(Sd_2)+ Uobs * Uobs * 1/Sd_2) every entries are scalars*/
- rResult[0] = rResult[0] -0.5 * (log(rResult[1])+ Uobs * Uobs /rResult[1]);
+ rResult[0] += -0.5 * (log(rResult[1])+ Uobs * Uobs /rResult[1]);
/* manual debug */
-
-
}
- UNPROTECT(11);
+ UNPROTECT(10);
return Result;
- }
-
- /* t=1;*/
-
-
- /* next step of the loop */
-
-
- /* for(i=0; i<p; i++){
- rMat21[i] = 0;
- for(j=0; j<p; j++){
- /* statevar <- expA %*% statevar: px1=pxp px1 */
- /* rMat21[i] = rMat21[i]+rExpA[i+j*p]*rStateVar[j];
- }
- rStateVar[i] = rMat21[i];
- rMat21[i]=0;/* statevar <- expA %*% statevar */
- /* }
-
-
-
- /* SigMatr <- expA %*% SigMatr %*% expAT + Qmatr: pxp = pxp pxp pxp */
- /* First We compute rMat22int <- expA %*% SigMatr : pxp = pxp pxp */
- /* for(i=0; i<p; i++){
- for(j=0; j<p; j++){
- rMat22int[i+j*p]=0;
- for(h=0; h<p; h++){
- rMat22int[i+j*p]=rMat22int[i+j*p]+rExpA[i+h*p]*
- rSigMatr[h+j*p];
- }
- }
- }
-
- /* print(rMat22int) */
- /* Second We compute rMat22est <- rMat22int %*% t(expA) + Qmatr: pxp = pxp pxp + pxp */
- /* for(i=0; i<p; i++){
- for(j=0; j<p; j++){
- rMat22est[i+j*p]=0;
- for(h=0; h<p; h++){
- rMat22est[i+j*p]=rMat22est[i+j*p]+rMat22int[i+h*p]*rExpA[j+h*p];
-
- }
- rSigMatr[i+j*p]=rMat22est[i+j*p];
-
- }
-
- }
-
-
-
-
- for(i=0; i<p*p; i++){
- rSigMatr[i]=rSigMatr[i]+rQmatr[i];
- }
-
-
- for(i=0; i<p*p; i++){
- rMat22est[i]=0;
- rMat22int[i]=0;
- }
-
- /*forecast*/
-
- /*Uobs <- y[t] - zc %*% statevar # 1 - 1Xp px1*/
- /* rMat22est[0]=0;
- for(i=0; i<p; i++){
- rMat22est[0]=rMat22est[0]+rZc[i]*rStateVar[i];
- }
- Uobs=rY[t]-rMat22est[0];
- rMat22est[0]=0;
- /* dum.zc <- zc %*% SigMatr # 1xp pxp */
-
-
-
- /*
-
-
-
- for(i=0; i<p; i++){
- rResult[1] = rResult[1]+rdum_zc[i]*rZc[i];
- }*/
-
-
-
- /* for(i=0; i<p; i++){
- rMat21[i] = 0;
- for(j=0; j<p; j++){
- rMat21[i] =rMat21[i]+rSigMatr[i+j*p]*rZc[j];
- }
- rdum_zc[i]=rMat21[i];
- rMat21[i]=0;
- }
-
- /*rResult[1]=rSigMatr[0];*/
- /* Sd_2 <- dum.zc %*% zcT # 1xp px1 */
- /* rResult[1] =0;
- for(i=0; i<p; i++){
- rResult[1] = rResult[1]+rdum_zc[i]*rZc[i];
- }
-
- /* correction */
- /* Kgain <- SigMatr %*% zcT %*% 1/Sd_2 # pxp px1*1 */
- /* for(i=0; i<p; i++){
- rMat21[i]=0;
- for(j=0; j<p; j++){
- rMat21[i]=rMat21[i]+rSigMatr[i+j*p]*rZc[j];
- }
- rKgain[i]=rMat21[i]/rResult[1];
- rMat21[i]=0;
- }
-
-
- /* statevar <- statevar+Kgain %*% Uobs # px1+px1 */
- /* for(i=0; i<p; i++){
- rStateVar[i] = rStateVar[i] + rKgain[i]*Uobs;
- }
-
-
- /* SigMatr <- SigMatr - Kgain %*% dum.zc # pxp-px1 1xp */
- /* for(i=0; i<p; i++){
- for(j=0; j<p; j++){
- rSigMatr[i+j*p]=rSigMatr[i+j*p]-rKgain[i]*rdum_zc[j];
- }
- }
- /*for(i=0; i<p*p; i++){
- rMat22int[i];
- }*/
- /*term_int = -0.5 * (log(Sd_2)+ Uobs * Uobs * 1/Sd_2) every entries are scalars*/
- /* rResult[0] = rResult[0] -0.5 * (log(rResult[1])+ Uobs * Uobs /rResult[1]);
-
-
-
-
- /* Unprotect steps */
- /* UNPROTECT(12);
- return(Result);
-
-
-
-
-
-
}
+
/*
SEXP sde_sim_euler(SEXP x0, SEXP t0, SEXP delta, SEXP N, SEXP M,
SEXP d,888 SEXP s, SEXP sx,
More information about the Yuima-commits
mailing list