[Yuima-commits] r336 - in pkg/yuima: . R

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

Author: iacus
Date: 2014-10-13 16:57:11 +0200 (Mon, 13 Oct 2014)
New Revision: 336

Modified:
   pkg/yuima/DESCRIPTION
   pkg/yuima/R/qmle.R
Log:
carma C

Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION	2014-09-29 12:25:28 UTC (rev 335)
+++ pkg/yuima/DESCRIPTION	2014-10-13 14:57:11 UTC (rev 336)
@@ -1,8 +1,8 @@
 Package: yuima
 Type: Package
 Title: The YUIMA Project package for SDEs
-Version: 1.0.36
-Date: 2014-09-29
+Version: 1.0.37
+Date: 2014-10-13
 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-09-29 12:25:28 UTC (rev 335)
+++ pkg/yuima/R/qmle.R	2014-10-13 14:57:11 UTC (rev 336)
@@ -1880,7 +1880,22 @@
   
   expAT<-t(expA)
   #SIGMA_err<-V_inf-expA%*%V_inf%*%t(expA)
-  SigMatr<-V_inf-expA%*%V_inf%*%expAT
+  
+  # input: V_inf, p, expA
+  # output: SigMatr (pre-alloc in R)
+  
+  SigMatr1 <- .Call("carma_tmp",  V_inf, as.integer(p), expA, PACKAGE="yuima")
+
+print( (expA %*% V_inf) %*% expAT)
+  SigMatr <- V_inf - expA %*% V_inf %*% expAT
+  cat("\nSigMatr\n")
+  print(SigMatr)
+  cat("\nAVAT\n")
+  print(expA %*% V_inf)
+  cat("\nSigMatr1\n")
+  print(SigMatr1)
+  stop("")
+  
   statevar<-matrix(rep(0, p),p,1)
   Qmatr<-SigMatr
   
@@ -1890,31 +1905,31 @@
   # SigMatr<-expA%*%V_inf%*%t(expA)+Qmatr
   
   #SigMatr<-Qmatr
-  SigMatr<-V_inf
+  SigMatr <- V_inf
   
-  zc<-matrix(bvector,1,p)
+  
+  zc <- matrix(bvector,1,p)
   loglstar <- 0
   loglstar1 <- 0
   
 #  zcT<-matrix(bvector,p,1)
-  zcT<-t(zc)
+  zcT <- t(zc)
+  ###  statevar, expA, times.obs, Qmatr,
   for(t in 1:times.obs){ 
     # prediction
-    statevar<-expA%*%statevar
-    SigMatr<-expA%*%SigMatr%*%expAT+Qmatr
+    statevar <- expA %*% statevar
+    SigMatr <- expA %*% SigMatr %*% expAT + Qmatr
     # forecast
-    Uobs<-y[t]-zc%*%statevar
-    dum.zc<-zc%*%SigMatr
-    sd_2<-dum.zc%*%zcT
-    # sd_2<-zc%*%SigMatr%*%zcT
-    Inv_sd_2<-1/sd_2
+    Uobs <- y[t] - zc %*% statevar
+    dum.zc <- zc %*% SigMatr
+    sd_2 <- dum.zc %*% zcT
+    Inv_sd_2 <- 1/sd_2
     #correction
-    Kgain<-SigMatr%*%zcT%*%Inv_sd_2    
-    statevar<-statevar+Kgain%*%Uobs
-    #SigMatr<-SigMatr-Kgain%*%zc%*%SigMatr
-    SigMatr<-SigMatr-Kgain%*%dum.zc
-    term_int<--0.5*(log(sd_2)+Uobs%*%Uobs%*%Inv_sd_2)
-    loglstar<-loglstar+term_int
+    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)
+    loglstar <- loglstar + term_int
   }
   return(list(loglstar=(loglstar-0.5*log(2*pi)*times.obs),s2hat=sd_2))
 }