[Yuima-commits] r603 - in pkg/yuima: . R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Apr 19 23:18:09 CEST 2017
Author: lorenzo
Date: 2017-04-19 23:18:08 +0200 (Wed, 19 Apr 2017)
New Revision: 603
Added:
pkg/yuima/src/Makevars
pkg/yuima/src/Makevars.win
pkg/yuima/src/pseudoLogLikCogarchIrregularGrid.cpp
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NAMESPACE
pkg/yuima/R/cogarchNoise.R
Log:
Added qmle on irregular grid for cogarch
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2017-04-12 11:19:28 UTC (rev 602)
+++ pkg/yuima/DESCRIPTION 2017-04-19 21:18:08 UTC (rev 603)
@@ -9,6 +9,7 @@
Description: Simulation and Inference for SDEs and Other Stochastic Processes.
License: GPL-2
URL: http://www.yuima-project.com
-LinkingTo: Rcpp
+LinkingTo: Rcpp, RcppArmadillo
+
Modified: pkg/yuima/NAMESPACE
===================================================================
--- pkg/yuima/NAMESPACE 2017-04-12 11:19:28 UTC (rev 602)
+++ pkg/yuima/NAMESPACE 2017-04-19 21:18:08 UTC (rev 603)
@@ -17,8 +17,9 @@
#08/07/2016
#exportPattern("^[[:alpha:]]+") # NEVER DO THIS AGAIN PLEASE!
+#cimport(RcppArmadillo)
importFrom(Rcpp, evalCpp)
-
+#import(RcppArmadillo)
# 03/07/2015
importFrom(stats, time)
importFrom(stats, ts)
@@ -198,8 +199,8 @@
export(lasso)
export(CPoint)
export(qmleR)
-export(qmleL)
-export(qmleLevy)
+export(qmleL)
+export(qmleLevy)
export(IC)
Modified: pkg/yuima/R/cogarchNoise.R
===================================================================
--- pkg/yuima/R/cogarchNoise.R 2017-04-12 11:19:28 UTC (rev 602)
+++ pkg/yuima/R/cogarchNoise.R 2017-04-19 21:18:08 UTC (rev 603)
@@ -59,7 +59,8 @@
b <- param[ar.name]
cost<- param[loc.par]
Data<-as.matrix(onezoo(observ)[,1])
- freq<-round(frequency(onezoo(observ)[,1]))
+ #freq<-round(frequency(onezoo(observ)[,1]))
+ freq<- diff(time(onezoo(observ)))
res<-auxcogarch.noise(cost,b,acoeff,mu,Data,freq,model at solve.variable)
return(res)
}
@@ -81,9 +82,11 @@
# Process_Y1 <- ExpY0
# Process_Y <- as.matrix(50.33)
var_V<-cost + sum(acoeff*Process_Y)
- delta <- 1/freq
+ # delta <- 1/freq
+ deltatot <- c(0,freq)
for(t in c(2:(length(Data)))){
# Y_t=e^{A\Delta t}Y_{t-\Delta t}+e^{A\left(\Delta t\right)}e\left(\Delta G_{t}\right)^{2}
+ delta <- deltatot[t]
Process_Y <- cbind(Process_Y, (expm(B*delta)%*%(Process_Y[,t-1]+e*squaredG[t])))
# Process_Y1 <- merge(Process_Y1, (expm(B*delta)%*%(Process_Y1[,t-1]+e*squaredG[t])))
#Process_Y <- cbind(Process_Y, (Process_Y[,t-1]+delta*B%*%Process_Y[,t-1]+e*squaredG[t]))
@@ -99,7 +102,10 @@
DumRet <- as.numeric(Data[1]+cumsum(DeltaG))
Cog <- cbind(DumRet,CogDum)
colnames(Cog) <- lab
- Cogarch <- setData(Cog, delta = delta)
+ Cog <- zoo(Cog,
+ order.by = cumsum(deltatot)) # Dataset on true irregular grid
+
+ Cogarch <- setData(Cog)
res<-list(incr.L=incr.L, Cogarch=Cogarch)
return(res)
}
Added: pkg/yuima/src/Makevars
===================================================================
--- pkg/yuima/src/Makevars (rev 0)
+++ pkg/yuima/src/Makevars 2017-04-19 21:18:08 UTC (rev 603)
@@ -0,0 +1 @@
+PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
Added: pkg/yuima/src/Makevars.win
===================================================================
--- pkg/yuima/src/Makevars.win (rev 0)
+++ pkg/yuima/src/Makevars.win 2017-04-19 21:18:08 UTC (rev 603)
@@ -0,0 +1,2 @@
+
+PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
Added: pkg/yuima/src/pseudoLogLikCogarchIrregularGrid.cpp
===================================================================
--- pkg/yuima/src/pseudoLogLikCogarchIrregularGrid.cpp (rev 0)
+++ pkg/yuima/src/pseudoLogLikCogarchIrregularGrid.cpp 2017-04-19 21:18:08 UTC (rev 603)
@@ -0,0 +1,81 @@
+#include <RcppArmadillo.h>
+
+using namespace Rcpp;
+
+/*
+Inputs
+int lengthObs
+arma::mat B, arma::mat Btilde, arma::mat InvBtilde
+double a0, double bq, double a1, double V, double PseudologLik
+arma::rowvec ta
+
+arma::colvec state, arma::colvec stateMean, arma::colvec e
+arma::vec DeltaG2, arma::vec Deltat
+
+
+double Irregular_PseudoLoglik_COG(int lengthObs, arma::mat B, arma::mat Btilde, arma::mat InvBtilde,
+ double a0, double bq, double a1, double V, double PseudologLik,arma::rowvec ta,
+ arma::colvec state, arma::colvec stateMean, arma::colvec e,
+ arma::vec DeltaG2, arma::vec Deltat) {
+ double res=0;
+ double VarDeltaG = 0;
+ arma::mat I = arma::eye<arma::mat>(B.n_rows, B.n_rows);
+ // for(i in c(1:(lengthObs))){
+ for(int i=0;i<lengthObs;i++){
+ arma::mat DeltatB1 = expmat(B * Deltat[i]);
+ arma::mat DeltatB2 = expmat(Btilde * Deltat[i]);
+ VarDeltaG = a0 * arma::as_scalar(Deltat[i]) * bq / (bq - a1); //+ arma::as_scalar(ta * DeltatB2 * InvBtilde * (I - DeltatB2) * (state - stateMean));
+ state <- (I + DeltaG2[i] / V * e * ta) * DeltatB1 * state+a0 * DeltaG2[i] / V * e;
+ if(VarDeltaG>0){
+ res =res - 0.5*(arma::as_scalar((DeltaG2[i] / VarDeltaG )) + log(VarDeltaG) + log(2 * 3.141593));
+ // res = res - 1/2 * (arma::as_scalar(DeltaG2[i] / VarDeltaG) + log(VarDeltaG) + log(2 * 3.141593));
+ }else{
+ res = res - 1000000.6543;
+ }
+ }
+ return PseudologLik;
+}
+*/
+
+// [[Rcpp::export]]
+double Irregular_PseudoLoglik_COG(int lengthObs, arma::mat B, arma::mat Btilde, arma::mat InvBtilde,
+ double a0, double bq, double a1, double V, double PseudologLik,arma::rowvec ta,
+ arma::colvec state, arma::colvec stateMean, arma::colvec e,
+ arma::vec DeltaG2, arma::vec Deltat) {
+ double res = 0;
+ double VarDeltaG = 0;
+ double penal = 100000.123;
+ arma::mat I = arma::eye<arma::mat>(B.n_rows, B.n_rows);
+// for(i in c(1:(lengthObs))){
+ for(int i=0;i<lengthObs;i++){
+ arma::mat DeltatB1 = expmat(B * Deltat[i]);
+ arma::mat DeltatB2 = expmat(Btilde * Deltat[i]);
+ VarDeltaG = a0 * Deltat[i] * bq / (bq - a1) + arma::as_scalar(ta * DeltatB2 * InvBtilde * (I - DeltatB2) * (state - stateMean));
+ state <- (I + DeltaG2[i] / V * e * ta) * DeltatB1 * state+a0 * DeltaG2[i] / V * e;
+ if(VarDeltaG>0){
+ res =res - 0.5*(arma::as_scalar((DeltaG2[i] / VarDeltaG )) + log(VarDeltaG) + log(2 * 3.141593));//
+ }else{
+ res = res+VarDeltaG-penal;
+ }
+}
+return res;
+}
+
+
+/*
+for(i in c(1:(lengthObs))){
+DeltatB1 <- expm(B*Deltat[i])
+DeltatB2 <- expm(Btilde*Deltat[i])
+VarDeltaG <- as.numeric(a0*Deltat[i]*bq/(bq-a1)+ta%*%DeltatB2%*%InvBtilde%*%(I-DeltatB2)%*%(state-stateMean))
+state <- (I+DeltaG2[i]/V*e%*%ta)%*%DeltatB1%*%state+a0*DeltaG2[i]/V*e
+V <- as.numeric(a0+ta%*% state)
+if(VarDeltaG>0){
+PseudologLik<--1/2*(DeltaG2[i]/VarDeltaG+log(VarDeltaG)+log(2*pi))+PseudologLik
+}else{
+PseudologLik<-PseudologLik-10^6
+}
+
+
+}
+*/
+
More information about the Yuima-commits
mailing list