[Yuima-commits] r510 - in pkg/yuima: . src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Nov 6 15:24:33 CET 2016
Author: iacus
Date: 2016-11-06 15:24:33 +0100 (Sun, 06 Nov 2016)
New Revision: 510
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/src/PseudoLoglikCOGARCH.c
pkg/yuima/src/RcppExports.cpp
pkg/yuima/src/carmafilter.c
pkg/yuima/src/cce_functions.c
pkg/yuima/src/qmlecpp.cpp
Log:
fixed line width
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2016-11-03 19:02:29 UTC (rev 509)
+++ pkg/yuima/DESCRIPTION 2016-11-06 14:24:33 UTC (rev 510)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project Package for SDEs
-Version: 1.3.5
+Version: 1.3.6
Depends: R(>= 2.10.0), methods, zoo, stats4, utils, expm, cubature, mvtnorm
Imports: Rcpp (>= 0.12.1)
Author: YUIMA Project Team
Modified: pkg/yuima/src/PseudoLoglikCOGARCH.c
===================================================================
--- pkg/yuima/src/PseudoLoglikCOGARCH.c 2016-11-03 19:02:29 UTC (rev 509)
+++ pkg/yuima/src/PseudoLoglikCOGARCH.c 2016-11-06 14:24:33 UTC (rev 510)
@@ -1,294 +1,294 @@
-/* The function computes the pseudo-loglikelihood of a COGARCH(p,q) model
- * See Iacus et al. 2015 for details
- */
-#include <Rmath.h>
-#include <R_ext/Boolean.h>
-#include <R_ext/Rdynload.h>
-#include <Rdefines.h>
-#include <Rinternals.h>
-#include <R_ext/Complex.h>
-
-#define max(a, b) (a > b ? a : b)
-#define min(a, b) (a < b ? a : b)
-
-/*SEXP myfun1(SEXP a, SEXP b);*/
-
-SEXP pseudoLoglik_COGARCH1(SEXP a0, SEXP bq, SEXP a1, SEXP stateMean, SEXP Q,
- SEXP DeltaG2, SEXP Deltat, SEXP DeltaB1,
- SEXP a, SEXP e,
- SEXP V, SEXP nObs,
- SEXP dummyMatr, SEXP dummyeB1);
-
-/*SEXP myfun1(SEXP a, SEXP b){
- double *ra, *rb, *rvai;
- double dummy1=0;
- int i, j;
- SEXP vai;
-
- PROTECT(a = AS_NUMERIC(a));
- ra = REAL(a);
-
- PROTECT(b = AS_NUMERIC(b));
- rb = REAL(b);
-
- PROTECT(vai = NEW_NUMERIC(1));
- rvai = REAL(vai);
-
- printf("\n Q %f.5", dummy1);
- for(i=0; i < 2; i++){
- dummy1 = 0;
- for(j=0; j < 2; j++){
- dummy1 += ra[i+j*2]*rb[j];
- printf("\n Q c: %d, %d %f.5", i, j, dummy1);
- }
- rvai[i]= dummy1;
-
- }
- UNPROTECT(3);
- return vai;
- }
-*/
-
-
-
-SEXP pseudoLoglik_COGARCH1(SEXP a0, SEXP bq, SEXP a1, SEXP stateMean, SEXP Q,
- SEXP DeltaG2, SEXP Deltat, SEXP DeltaB1,
- SEXP a, SEXP e,
- SEXP V, SEXP nObs,
- SEXP dummyMatr, SEXP dummyeB1){
-
- /* Declare Integer Variable */
- int numb_Obs, q, t, i, j;
- double *ra0, *rbq, *ra1, *rstateMean, *rstate;
- double *rDeltaG2, *rDeltat, *rDeltaB1;
- double *ra, *re, *rV, rVarDeltaG=0;
- double *rdummyMatr;
- /*rPseudologLik=0,*/
- double *rdummyeB1;
- double dummy1=0;
- double dummy2=0;
- double *res;
- double *rstatedum;
- SEXP PseudoLoglik;
- SEXP state;
- SEXP statedum;
-
-
-
-
- /* Protect Objects */
-
- PROTECT(a0 = AS_NUMERIC(a0));
- ra0 = REAL(a0);
-
- /*1*/
-
- PROTECT(PseudoLoglik = NEW_NUMERIC(1));
- res=REAL(PseudoLoglik);
-
-
-
- /*2*/
-
- PROTECT(bq = AS_NUMERIC(bq));
- rbq = REAL(bq);
-
- /*3*/
-
- PROTECT(a1 = AS_NUMERIC(a1));
- ra1 = REAL(a1);
-
- /*4*/
-
- PROTECT(stateMean = AS_NUMERIC(stateMean));
- rstateMean = REAL(stateMean);
-
- /*5*/
-
-
-
- PROTECT(DeltaG2 = AS_NUMERIC(DeltaG2));
- rDeltaG2 = REAL(DeltaG2);
-
- /*6*/
-
- PROTECT(Deltat = AS_NUMERIC(Deltat));
- rDeltat = REAL(Deltat);
-
- /*7*/
-
- PROTECT(DeltaB1 = AS_NUMERIC(DeltaB1));
- rDeltaB1 = REAL(DeltaB1);
-
- /*8*/
-
- PROTECT(a = AS_NUMERIC(a));
- ra = REAL(a);
-
- /*9*/
-
-
- PROTECT(e = AS_NUMERIC(e));
- re = REAL(e);
-
- /*10*/
-
-
- PROTECT(V = AS_NUMERIC(V));
- rV = REAL(V);
-
-
-
-
- /*11*/
-
-
- PROTECT(dummyMatr = AS_NUMERIC(dummyMatr));
- rdummyMatr = REAL(dummyMatr);
-
- /*12*/
-
- PROTECT(dummyeB1 = AS_NUMERIC(dummyeB1));
- rdummyeB1 = REAL(dummyeB1);
-
- /*13*/
-
-
- /* Declare dimensions for the state variables and observations */
- numb_Obs = *INTEGER(nObs);
- q = *INTEGER(Q);
-
- PROTECT(state = NEW_NUMERIC(q));
- rstate = REAL(state);
-
- PROTECT(statedum = NEW_NUMERIC(q));
- rstatedum = REAL(statedum);
-
- for(i=0; i<q; i++){
- rstate[i]=rstateMean[i];
- }
-
- /*printf("\n Q c: %d, %d ", q, numb_Obs);
- for(i=0; i<q; i++){
- printf("\n RSTATE: %.5f %.5f %.5f", rstate[i], rstateMean[i], rdummyMatr[i]);
- printf("\n RMAtr: %.5f %.5f" , rDeltaB1[i], rdummyeB1[i]);
- }*/
-
- *res=0;
- /*fd = fopen("dueinteri.txt", "r");
- printf("\n %p %p", &rstate, &rstateMean);*/
- for(t=0; t<numb_Obs; t++){
- /* VarDeltaG <- as.numeric(a0*Deltat*bq/(bq-a1)+dummyMatr%*%(state-stateMean)) */
- dummy1 = 0;
- for(i=0; i<q; i++){
- dummy1 += rdummyMatr[i]*(rstate[i]-rstateMean[i]);
- /*printf("\n %d %.10f %.10f %.10f %.10f", i, dummy1, rdummyMatr[i], rstate[i], rstateMean[i]);*/
- }
-
- /*printf("\n dummy1 c: %.10f", dummy1);*/
-
- rVarDeltaG = ra0[0]*rDeltat[0]*rbq[0]/(rbq[0]-ra1[0])+dummy1;
-
- /* state <- DeltatB1%*%state+DeltaG2[i]/V*dummyeB1%*%state+a0*DeltaG2[i]/V*e */
-
- for(i=0; i<q; i++){
- dummy1 = 0;
- dummy2 = 0;
- for(j=0; j<q; j++){
- dummy1 += rDeltaB1[i+j*q]*rstate[j];
- dummy2 += rdummyeB1[i+j*q]*rstate[j];
-
- }
- rstatedum[i]= dummy1+rDeltaG2[t]/rV[0]*dummy2+ra0[0]*rDeltaG2[t]/rV[0]*re[i];
- /*rstate[i]= dummy1+dummy2;
- printf("\n d1 %.10f d2 %.10f", dummy1, dummy2);*/
- }
-
- for(i=0; i<q; i++){
- rstate[i]=rstatedum[i];
- }
- /* V <- as.numeric(a0+ta%*% state) */
- rV[0] = 0;
- for(i=0; i<q; i++){
- rV[0] += ra[i]*rstate[i];
- }
- rV[0] = rV[0]+ra0[0];
-
- /* PseudologLik<- -1/2*(DeltaG2[i]/VarDeltaG+log(VarDeltaG)+log(2*pi)) */
-
- *res += 0.5*(-rDeltaG2[t]/rVarDeltaG-log(rVarDeltaG)-log(2.*3.14159265));
- /*fscanf(fd, "%lf", &res);
- printf("\n c: %.10f - %.5f %.5f - %.5f", rVarDeltaG, rstate[0], rstate[1], rV[0]);*/
-
- }
-
- /*fclose(fd);*/
-
-
- UNPROTECT(15);
- return PseudoLoglik;
-
- }
-
-
-/*SEXP pseudoLoglik_COGARCH(SEXP a0, SEXP bq, SEXP a1, SEXP stateMean, SEXP Q,
- SEXP state, SEXP DeltaG2, SEXP Deltat, SEXP DeltatB,
- SEXP B, SEXP a,SEXP e, SEXP Btilde,
- SEXP InvBtilde, SEXP V, SEXP I, SEXP VarDeltaG,
- SEXP PseudologLik, SEXP nObs, SEXP expr, SEXP rho){
-
- Declare Integer Variable
-
-
-
- int numb_Obs, q, t;
- SEXP ans, A;
- double rA, rans;
-
- if(!isNewList(DeltatB)) error("'DeltatB' must be a list");
-
- if(!isEnvironment(rho)) error("`rho' must be an environment");
- Protect Objects
-
- numb_Obs = *INTEGER(nObs);
-
- q = *INTEGER(Q);
-
- PROTECT(B = allocMatrix(REALSXP, q, q));
- PROTECT(DeltatB);
-
- ans = PROTECT(allocVector(VECSXP, numb_Obs));
-
- PROTECT(A = allocMatrix(REALSXP, q, q));
- rA = *REAL(A);
- rPseudologLik[0]=0;
- for(t=0; t<numb_Obs; t++){
-
- VarDeltaG <- as.numeric(a0*Deltat[i]*bq/(bq-a1)+t(a)%*
- %expm(Btilde*Deltat[i]) %*%InvBtilde%*% (I-expm(-Btilde*Deltat[i]))
- %*%(state-stateMean))
-
-
-
- state <- (I+DeltaG2[i]/V*e%*%t(a))%*%expm(B*Deltat[i])
- %*%state+a0*DeltaG2[i]/V*e
-
-
- V <- as.numeric(a0+t(a)%*% state)
-
-
- PseudologLik<--1/2*(DeltaG2[i]/VarDeltaG+
- log(VarDeltaG)+log(2*pi))
-
- defineVar(install("x"), VECTOR_ELT(DeltatB,t), rho);
- SET_VECTOR_ELT(ans, t, eval(expr, rho));
-
- }
-
-
- UNPROTECT(4);
- return PseudologLik;
- return ans;
-
- }
- */
+/* The function computes the pseudo-loglikelihood of a COGARCH(p,q) model
+ * See Iacus et al. 2015 for details
+ */
+#include <Rmath.h>
+#include <R_ext/Boolean.h>
+#include <R_ext/Rdynload.h>
+#include <Rdefines.h>
+#include <Rinternals.h>
+#include <R_ext/Complex.h>
+
+#define max(a, b) (a > b ? a : b)
+#define min(a, b) (a < b ? a : b)
+
+/*SEXP myfun1(SEXP a, SEXP b);*/
+
+SEXP pseudoLoglik_COGARCH1(SEXP a0, SEXP bq, SEXP a1, SEXP stateMean, SEXP Q,
+ SEXP DeltaG2, SEXP Deltat, SEXP DeltaB1,
+ SEXP a, SEXP e,
+ SEXP V, SEXP nObs,
+ SEXP dummyMatr, SEXP dummyeB1);
+
+/*SEXP myfun1(SEXP a, SEXP b){
+ double *ra, *rb, *rvai;
+ double dummy1=0;
+ int i, j;
+ SEXP vai;
+
+ PROTECT(a = AS_NUMERIC(a));
+ ra = REAL(a);
+
+ PROTECT(b = AS_NUMERIC(b));
+ rb = REAL(b);
+
+ PROTECT(vai = NEW_NUMERIC(1));
+ rvai = REAL(vai);
+
+ printf("\n Q %f.5", dummy1);
+ for(i=0; i < 2; i++){
+ dummy1 = 0;
+ for(j=0; j < 2; j++){
+ dummy1 += ra[i+j*2]*rb[j];
+ printf("\n Q c: %d, %d %f.5", i, j, dummy1);
+ }
+ rvai[i]= dummy1;
+
+ }
+ UNPROTECT(3);
+ return vai;
+ }
+*/
+
+
+
+SEXP pseudoLoglik_COGARCH1(SEXP a0, SEXP bq, SEXP a1, SEXP stateMean, SEXP Q,
+ SEXP DeltaG2, SEXP Deltat, SEXP DeltaB1,
+ SEXP a, SEXP e,
+ SEXP V, SEXP nObs,
+ SEXP dummyMatr, SEXP dummyeB1){
+
+ /* Declare Integer Variable */
+ int numb_Obs, q, t, i, j;
+ double *ra0, *rbq, *ra1, *rstateMean, *rstate;
+ double *rDeltaG2, *rDeltat, *rDeltaB1;
+ double *ra, *re, *rV, rVarDeltaG=0;
+ double *rdummyMatr;
+ /*rPseudologLik=0,*/
+ double *rdummyeB1;
+ double dummy1=0;
+ double dummy2=0;
+ double *res;
+ double *rstatedum;
+ SEXP PseudoLoglik;
+ SEXP state;
+ SEXP statedum;
+
+
+
+
+ /* Protect Objects */
+
+ PROTECT(a0 = AS_NUMERIC(a0));
+ ra0 = REAL(a0);
+
+ /*1*/
+
+ PROTECT(PseudoLoglik = NEW_NUMERIC(1));
+ res=REAL(PseudoLoglik);
+
+
+
+ /*2*/
+
+ PROTECT(bq = AS_NUMERIC(bq));
+ rbq = REAL(bq);
+
+ /*3*/
+
+ PROTECT(a1 = AS_NUMERIC(a1));
+ ra1 = REAL(a1);
+
+ /*4*/
+
+ PROTECT(stateMean = AS_NUMERIC(stateMean));
+ rstateMean = REAL(stateMean);
+
+ /*5*/
+
+
+
+ PROTECT(DeltaG2 = AS_NUMERIC(DeltaG2));
+ rDeltaG2 = REAL(DeltaG2);
+
+ /*6*/
+
+ PROTECT(Deltat = AS_NUMERIC(Deltat));
+ rDeltat = REAL(Deltat);
+
+ /*7*/
+
+ PROTECT(DeltaB1 = AS_NUMERIC(DeltaB1));
+ rDeltaB1 = REAL(DeltaB1);
+
+ /*8*/
+
+ PROTECT(a = AS_NUMERIC(a));
+ ra = REAL(a);
+
+ /*9*/
+
+
+ PROTECT(e = AS_NUMERIC(e));
+ re = REAL(e);
+
+ /*10*/
+
+
+ PROTECT(V = AS_NUMERIC(V));
+ rV = REAL(V);
+
+
+
+
+ /*11*/
+
+
+ PROTECT(dummyMatr = AS_NUMERIC(dummyMatr));
+ rdummyMatr = REAL(dummyMatr);
+
+ /*12*/
+
+ PROTECT(dummyeB1 = AS_NUMERIC(dummyeB1));
+ rdummyeB1 = REAL(dummyeB1);
+
+ /*13*/
+
+
+ /* Declare dimensions for the state variables and observations */
+ numb_Obs = *INTEGER(nObs);
+ q = *INTEGER(Q);
+
+ PROTECT(state = NEW_NUMERIC(q));
+ rstate = REAL(state);
+
+ PROTECT(statedum = NEW_NUMERIC(q));
+ rstatedum = REAL(statedum);
+
+ for(i=0; i<q; i++){
+ rstate[i]=rstateMean[i];
+ }
+
+ /*printf("\n Q c: %d, %d ", q, numb_Obs);
+ for(i=0; i<q; i++){
+ printf("\n RSTATE: %.5f %.5f %.5f", rstate[i], rstateMean[i], rdummyMatr[i]);
+ printf("\n RMAtr: %.5f %.5f" , rDeltaB1[i], rdummyeB1[i]);
+ }*/
+
+ *res=0;
+ /*fd = fopen("dueinteri.txt", "r");
+ printf("\n %p %p", &rstate, &rstateMean);*/
+ for(t=0; t<numb_Obs; t++){
+ /* VarDeltaG <- as.numeric(a0*Deltat*bq/(bq-a1)+dummyMatr%*%(state-stateMean)) */
+ dummy1 = 0;
+ for(i=0; i<q; i++){
+ dummy1 += rdummyMatr[i]*(rstate[i]-rstateMean[i]);
+ /*printf("\n %d %.10f %.10f %.10f %.10f", i, dummy1, rdummyMatr[i], rstate[i], rstateMean[i]);*/
+ }
+
+ /*printf("\n dummy1 c: %.10f", dummy1);*/
+
+ rVarDeltaG = ra0[0]*rDeltat[0]*rbq[0]/(rbq[0]-ra1[0])+dummy1;
+
+ /* state <- DeltatB1%*%state+DeltaG2[i]/V*dummyeB1%*%state+a0*DeltaG2[i]/V*e */
+
+ for(i=0; i<q; i++){
+ dummy1 = 0;
+ dummy2 = 0;
+ for(j=0; j<q; j++){
+ dummy1 += rDeltaB1[i+j*q]*rstate[j];
+ dummy2 += rdummyeB1[i+j*q]*rstate[j];
+
+ }
+ rstatedum[i]= dummy1+rDeltaG2[t]/rV[0]*dummy2+ra0[0]*rDeltaG2[t]/rV[0]*re[i];
+ /*rstate[i]= dummy1+dummy2;
+ printf("\n d1 %.10f d2 %.10f", dummy1, dummy2);*/
+ }
+
+ for(i=0; i<q; i++){
+ rstate[i]=rstatedum[i];
+ }
+ /* V <- as.numeric(a0+ta%*% state) */
+ rV[0] = 0;
+ for(i=0; i<q; i++){
+ rV[0] += ra[i]*rstate[i];
+ }
+ rV[0] = rV[0]+ra0[0];
+
+ /* PseudologLik<- -1/2*(DeltaG2[i]/VarDeltaG+log(VarDeltaG)+log(2*pi)) */
+
+ *res += 0.5*(-rDeltaG2[t]/rVarDeltaG-log(rVarDeltaG)-log(2.*3.14159265));
+ /*fscanf(fd, "%lf", &res);
+ printf("\n c: %.10f - %.5f %.5f - %.5f", rVarDeltaG, rstate[0], rstate[1], rV[0]);*/
+
+ }
+
+ /*fclose(fd);*/
+
+
+ UNPROTECT(15);
+ return PseudoLoglik;
+
+ }
+
+
+/*SEXP pseudoLoglik_COGARCH(SEXP a0, SEXP bq, SEXP a1, SEXP stateMean, SEXP Q,
+ SEXP state, SEXP DeltaG2, SEXP Deltat, SEXP DeltatB,
+ SEXP B, SEXP a,SEXP e, SEXP Btilde,
+ SEXP InvBtilde, SEXP V, SEXP I, SEXP VarDeltaG,
+ SEXP PseudologLik, SEXP nObs, SEXP expr, SEXP rho){
+
+ Declare Integer Variable
+
+
+
+ int numb_Obs, q, t;
+ SEXP ans, A;
+ double rA, rans;
+
+ if(!isNewList(DeltatB)) error("'DeltatB' must be a list");
+
+ if(!isEnvironment(rho)) error("`rho' must be an environment");
+ Protect Objects
+
+ numb_Obs = *INTEGER(nObs);
+
+ q = *INTEGER(Q);
+
+ PROTECT(B = allocMatrix(REALSXP, q, q));
+ PROTECT(DeltatB);
+
+ ans = PROTECT(allocVector(VECSXP, numb_Obs));
+
+ PROTECT(A = allocMatrix(REALSXP, q, q));
+ rA = *REAL(A);
+ rPseudologLik[0]=0;
+ for(t=0; t<numb_Obs; t++){
+
+ VarDeltaG <- as.numeric(a0*Deltat[i]*bq/(bq-a1)+t(a)%*
+ %expm(Btilde*Deltat[i]) %*%InvBtilde%*% (I-expm(-Btilde*Deltat[i]))
+ %*%(state-stateMean))
+
+
+
+ state <- (I+DeltaG2[i]/V*e%*%t(a))%*%expm(B*Deltat[i])
+ %*%state+a0*DeltaG2[i]/V*e
+
+
+ V <- as.numeric(a0+t(a)%*% state)
+
+
+ PseudologLik<--1/2*(DeltaG2[i]/VarDeltaG+
+ log(VarDeltaG)+log(2*pi))
+
+ defineVar(install("x"), VECTOR_ELT(DeltatB,t), rho);
+ SET_VECTOR_ELT(ans, t, eval(expr, rho));
+
+ }
+
+
+ UNPROTECT(4);
+ return PseudologLik;
+ return ans;
+
+ }
+ */
Modified: pkg/yuima/src/RcppExports.cpp
===================================================================
--- pkg/yuima/src/RcppExports.cpp 2016-11-03 19:02:29 UTC (rev 509)
+++ pkg/yuima/src/RcppExports.cpp 2016-11-06 14:24:33 UTC (rev 510)
@@ -1,67 +1,67 @@
-// This file was generated by Rcpp::compileAttributes
-// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
-
-#include <Rcpp.h>
-
-using namespace Rcpp;
-
-// detcpp
-double detcpp(NumericMatrix A);
-RcppExport SEXP yuima_detcpp(SEXP ASEXP) {
-BEGIN_RCPP
- Rcpp::RObject __result;
- Rcpp::RNGScope __rngScope;
- Rcpp::traits::input_parameter< NumericMatrix >::type A(ASEXP);
- __result = Rcpp::wrap(detcpp(A));
- return __result;
-END_RCPP
-}
-// Smake
-NumericMatrix Smake(NumericVector b, int d);
-RcppExport SEXP yuima_Smake(SEXP bSEXP, SEXP dSEXP) {
-BEGIN_RCPP
- Rcpp::RObject __result;
- Rcpp::RNGScope __rngScope;
- Rcpp::traits::input_parameter< NumericVector >::type b(bSEXP);
- Rcpp::traits::input_parameter< int >::type d(dSEXP);
- __result = Rcpp::wrap(Smake(b, d));
- return __result;
-END_RCPP
-}
-// solvecpp
-NumericMatrix solvecpp(NumericMatrix A);
-RcppExport SEXP yuima_solvecpp(SEXP ASEXP) {
-BEGIN_RCPP
- Rcpp::RObject __result;
- Rcpp::RNGScope __rngScope;
- Rcpp::traits::input_parameter< NumericMatrix >::type A(ASEXP);
- __result = Rcpp::wrap(solvecpp(A));
- return __result;
-END_RCPP
-}
-// trace
-double sub_f(NumericMatrix S, NumericVector b);
-RcppExport SEXP yuima_sub_f(SEXP SSEXP, SEXP bSEXP) {
-BEGIN_RCPP
- Rcpp::RObject __result;
- Rcpp::RNGScope __rngScope;
- Rcpp::traits::input_parameter< NumericMatrix >::type S(SSEXP);
- Rcpp::traits::input_parameter< NumericVector >::type b(bSEXP);
- __result = Rcpp::wrap(sub_f(S, b));
- return __result;
-END_RCPP
-}
-// likndim
-double likndim(NumericMatrix dx, NumericMatrix b, NumericMatrix A, double h);
-RcppExport SEXP yuima_likndim(SEXP dxSEXP, SEXP bSEXP, SEXP ASEXP, SEXP hSEXP) {
-BEGIN_RCPP
- Rcpp::RObject __result;
- Rcpp::RNGScope __rngScope;
- Rcpp::traits::input_parameter< NumericMatrix >::type dx(dxSEXP);
- Rcpp::traits::input_parameter< NumericMatrix >::type b(bSEXP);
- Rcpp::traits::input_parameter< NumericMatrix >::type A(ASEXP);
- Rcpp::traits::input_parameter< double >::type h(hSEXP);
- __result = Rcpp::wrap(likndim(dx, b, A, h));
- return __result;
-END_RCPP
-}
+// This file was generated by Rcpp::compileAttributes
+// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+#include <Rcpp.h>
+
+using namespace Rcpp;
+
+// detcpp
+double detcpp(NumericMatrix A);
+RcppExport SEXP yuima_detcpp(SEXP ASEXP) {
+BEGIN_RCPP
+ Rcpp::RObject __result;
+ Rcpp::RNGScope __rngScope;
+ Rcpp::traits::input_parameter< NumericMatrix >::type A(ASEXP);
+ __result = Rcpp::wrap(detcpp(A));
+ return __result;
+END_RCPP
+}
+// Smake
+NumericMatrix Smake(NumericVector b, int d);
+RcppExport SEXP yuima_Smake(SEXP bSEXP, SEXP dSEXP) {
+BEGIN_RCPP
+ Rcpp::RObject __result;
+ Rcpp::RNGScope __rngScope;
+ Rcpp::traits::input_parameter< NumericVector >::type b(bSEXP);
+ Rcpp::traits::input_parameter< int >::type d(dSEXP);
+ __result = Rcpp::wrap(Smake(b, d));
+ return __result;
+END_RCPP
+}
+// solvecpp
+NumericMatrix solvecpp(NumericMatrix A);
+RcppExport SEXP yuima_solvecpp(SEXP ASEXP) {
+BEGIN_RCPP
+ Rcpp::RObject __result;
+ Rcpp::RNGScope __rngScope;
+ Rcpp::traits::input_parameter< NumericMatrix >::type A(ASEXP);
+ __result = Rcpp::wrap(solvecpp(A));
+ return __result;
+END_RCPP
+}
+// trace
+double sub_f(NumericMatrix S, NumericVector b);
+RcppExport SEXP yuima_sub_f(SEXP SSEXP, SEXP bSEXP) {
+BEGIN_RCPP
+ Rcpp::RObject __result;
+ Rcpp::RNGScope __rngScope;
+ Rcpp::traits::input_parameter< NumericMatrix >::type S(SSEXP);
+ Rcpp::traits::input_parameter< NumericVector >::type b(bSEXP);
+ __result = Rcpp::wrap(sub_f(S, b));
+ return __result;
+END_RCPP
+}
+// likndim
+double likndim(NumericMatrix dx, NumericMatrix b, NumericMatrix A, double h);
+RcppExport SEXP yuima_likndim(SEXP dxSEXP, SEXP bSEXP, SEXP ASEXP, SEXP hSEXP) {
+BEGIN_RCPP
+ Rcpp::RObject __result;
+ Rcpp::RNGScope __rngScope;
+ Rcpp::traits::input_parameter< NumericMatrix >::type dx(dxSEXP);
+ Rcpp::traits::input_parameter< NumericMatrix >::type b(bSEXP);
+ Rcpp::traits::input_parameter< NumericMatrix >::type A(ASEXP);
+ Rcpp::traits::input_parameter< double >::type h(hSEXP);
+ __result = Rcpp::wrap(likndim(dx, b, A, h));
+ return __result;
+END_RCPP
+}
Modified: pkg/yuima/src/carmafilter.c
===================================================================
--- pkg/yuima/src/carmafilter.c 2016-11-03 19:02:29 UTC (rev 509)
+++ pkg/yuima/src/carmafilter.c 2016-11-06 14:24:33 UTC (rev 510)
@@ -1,380 +1,380 @@
-/*
- * R : A Computer Language for Statistical Data Analysis
- * Code of this package: Copyright (C) 2006 S. M. Iacus
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation; either version 2 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
- *
- *
- * Exports
- * sde_sim_xxx(...)
- *
- * to be called as .C(.) in ../R/sde.sim.xxx.R
- * where xxx is one among "euler", "milstein", "milstei2", "KPS"
- */
-
-
-
-#include <R.h>
-#include <Rmath.h>
-#include <R_ext/Boolean.h>
-#include <R_ext/Rdynload.h>
-#include <Rdefines.h>
-#include <Rinternals.h>
-#include <R_ext/Complex.h>
-
-#define max(a, b) (a > b ? a : b)
-#define min(a, b) (a < b ? a : b)
-
-
-
-
-
-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 Result,SEXP Kgain,
- SEXP dum_zc, SEXP Mat22int);
-
-
-
-SEXP carma_tmp(SEXP V, SEXP P, SEXP A){
-
- int p;
- int i, j, h;
- double *rV, *rA, *rB, *rC, *rSigma;
- SEXP B, C, Sigma;
-
- if(!isInteger(P)) error("`P' must be integer");
- if(!isNumeric(V)) error("`V' must be numeric");
- if(!isNumeric(A)) error("`A' must be numeric");
-
-
- PROTECT(V = AS_NUMERIC(V));
- rV = REAL(V);
- PROTECT(A = AS_NUMERIC(A));
- rA = REAL(A);
-
- p = *INTEGER(P);
-
-
- PROTECT(B = allocMatrix(REALSXP, p, p));
- rB = REAL(B);
-
- PROTECT(C = allocMatrix(REALSXP, p, p));
- rC = REAL(C);
-
-
- PROTECT(Sigma = allocMatrix(REALSXP, p, p));
- rSigma = REAL(Sigma);
-
- /* B = A %*% V */
- for(i=0; i<p; i++){
- for(j=0; j<p; j++){
- rB[i+j*p] = 0;
- for(h=0; h<p; h++){
- rB[i+j*p] = rB[i+j*p] + rA[i+h*p] * rV[h+j*p];
- }
- }
- /*}
- C = B %*% A^T
- for(i=0; i<p; i++){*/
- for(j=0; j<p; j++){
- rC[i+j*p] = 0;
- for(h=0; h<p; h++){
- rC[i+j*p] = rC[i+j*p] + rB[i+h*p] * rA[j+h*p];
- }
- rSigma[i+j*p] = rV[i+j*p] - rC[i+j*p];
- }
- }
- /*for(i=0; i<p*p; i++){
-
- }*/
- UNPROTECT(5);
- return Sigma;
-}
-
-SEXP Cycle_Carma(SEXP Y, SEXP StateVar, SEXP ExpA, SEXP Times_Obs, SEXP P,
- 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, t;
-
-
- /* Declare pointer */
-
- double *rY, *rStateVar, *rExpA, *rQmatr, *rSigMatr, *rZc;
- double *rKgain, *rdum_zc, *rMat22int, *rResult;
- double Uobs=0;
- double dummy=0;
- double dummy1=0;
- /*FILE *fd;*/
- /* Declare SEXP */
-
-
-
- /* 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");
-
-
- /* Protect Objects */
-
- PROTECT(Y = AS_NUMERIC(Y));
- rY = REAL(Y);
-
- PROTECT(StateVar = AS_NUMERIC(StateVar));
- rStateVar = REAL(StateVar);
-
- PROTECT(ExpA = AS_NUMERIC(ExpA));
- rExpA = REAL(ExpA);
-
- PROTECT(Qmatr = AS_NUMERIC(Qmatr));
- rQmatr = REAL(Qmatr);
-
- PROTECT(SigMatr = AS_NUMERIC(SigMatr));
- rSigMatr = REAL(SigMatr);
-
- 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);*/
-
-
- /* Declare dimensions for the state variables and observations */
- times_obs = *INTEGER(Times_Obs);
- p = *INTEGER(P);
-
-
- /* Main Code
- Dimension of Inputs:
- Y = Vector p dimension;
- StateVar = matrix p x 1;
- ExpA = matrix p x p;
- Times_Obs = integer;
- P = integer;
- Qmatr = matrix p x p;
- SigMatr = matrix p x p;
- Zc = matrix 1 x p;
- Logstar = scalar.
-
- for(t in 1:rtime_obs){*/
-
- rResult[0]=0;
-
- for(t=0; t<times_obs; t++){
- /*t=0;*/
- /* prediction */
- for(i=0; i<p; i++){
- dummy = 0;
- for(j=0; j<p; j++){
- /* statevar <- expA %*% statevar: px1=pxp px1 */
- dummy += rExpA[i+j*p]*rStateVar[j];
- }
- rStateVar[i] = dummy;
- /*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];
- }
- }
- }
- /* 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++){
- rSigMatr[i+j*p]=0;
- for(h=0; h<p; h++){
- rSigMatr[i+j*p] += rMat22int[i+h*p]*rExpA[j+h*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];
- }
- 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;
- /* Sd_2 <- dum.zc %*% zcT # 1xp px1 */
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 510
More information about the Yuima-commits
mailing list