[Yuima-commits] r365 - in pkg/yuima: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Mar 9 10:56:52 CET 2015
Author: lorenzo
Date: 2015-03-09 10:56:52 +0100 (Mon, 09 Mar 2015)
New Revision: 365
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/R/MM.COGARCH.R
pkg/yuima/R/cogarchNoise.R
pkg/yuima/R/simulate.R
Log:
Modify gmm cogarch
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2015-03-01 21:23:47 UTC (rev 364)
+++ pkg/yuima/DESCRIPTION 2015-03-09 09:56:52 UTC (rev 365)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package for SDEs
-Version: 1.0.57
-Date: 2015-02-19
+Version: 1.0.58
+Date: 2015-03-09
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/MM.COGARCH.R
===================================================================
--- pkg/yuima/R/MM.COGARCH.R 2015-03-01 21:23:47 UTC (rev 364)
+++ pkg/yuima/R/MM.COGARCH.R 2015-03-09 09:56:52 UTC (rev 365)
@@ -186,7 +186,7 @@
assign("d", d, envir=env)
typeacf <- "correlation"
assign("typeacf", typeacf, envir=env)
- CovQuad <- log(abs(acf(G_i^2,plot=FALSE,lag.max=min(d,env$lag),type=typeacf)$acf[-1]))
+ CovQuad <- acf(G_i^2,plot=FALSE,lag.max=min(d,env$lag),type=typeacf)$acf[-1]
#CovQuad <-log(abs(yuima.acf(data=G_i^2,lag.max=min(d,env$lag))))
assign("CovQuad", CovQuad, envir=env)
@@ -247,7 +247,7 @@
avect<-out$par[ma.name]
a1<-avect[1]
- out$par[loc.par]<-(bq-a1)*mu_G2/(bq*r)
+ #out$par[loc.par]<-(bq-a1)*mu_G2/(bq*r)
# Determine the Variance-Covariance Matrix
dimOutp<-length(out$par)-2
@@ -328,6 +328,7 @@
cost<-param[cost]
for(i in c(1:length(h))){
MomentCog <- MM_Cogarch(p = env$p, q = env$q, acoeff=param[a],
+ cost=cost,
b=param[b], r = r, h = h[i],
type = typeacf, m2=mu_G2, var=var_G2)
@@ -337,11 +338,11 @@
theo_mu_G2 <- MomentCog$meanG2
param[cost]<-MomentCog$cost
}
- res <- sum((c(log(abs(TheoCovQuad)))-c(CovQuad))^2)
+ res <- sum((c(theo_mu_G2,TheoCovQuad)-c(mu_G2,CovQuad))^2)
return(res)
}
-MM_Cogarch <- function(p, q, acoeff, b, r, h, type, m2, var){
+MM_Cogarch <- function(p, q, acoeff,cost, b, r, h, type, m2, var){
# The code developed here is based on the acf for squared returns derived in the
# Chaadra phd Thesis
a <- e <- matrix(0,nrow=q,ncol=1)
@@ -362,9 +363,9 @@
mu<-1 # we assume variance of the underlying L\'evy is one
meanL1<-mu
- cost<-(bq-mu*a1)*m2/(bq*r)
-
- B_tilde <- MatrixA(b[c(q:1)])+mu*e%*%t(a)
+ # cost<-(bq-mu*a1)*m2/(bq*r)
+ B<- MatrixA(b[c(q:1)])
+ B_tilde <- B+mu*e%*%t(a)
meanG2 <- cost*bq*r/(bq-mu*a1)*meanL1
Inf_eps <- IdentM <- diag(q)
if(q==1){
@@ -386,13 +387,19 @@
}
term <- expm(B_tilde*h)
- invB <- solve(B_tilde) # In this case we can use the analytical form for Companion Matrix???
+ #invB <- solve(B_tilde) # In this case we can use the analytical form for Companion Matrix???
+ # We invert the B_tilde using the Inverse of companion matrix
+ if(q>1){
+ invB<-rbind(c(-B_tilde[q,-1],1)/B_tilde[q,1],cbind(diag(q-1),matrix(0,q-1,1)))
+ }else{invB<-1/B_tilde}
term1 <- invB%*%(IdentM-expm(-B_tilde*r))
term2 <- (expm(B_tilde*r)-IdentM)
P0_overRho <- 2*mu^2*(3*invB%*%(invB%*%term2-r*IdentM)-IdentM)%*%Inf_eps
Q0_overRho <- 6*mu*((r*IdentM-invB%*%term2)%*%Inf_eps
-invB%*%(invB%*%term2-r*IdentM)%*%Inf_eps%*%t(B_tilde))%*%e
+# Q0_overRho <- 6*mu*((r*IdentM-invB%*%term2)%*%Inf_eps
+# -invB%*%(invB%*%term2-r*IdentM)%*%Inf_eps%*%t(B))%*%e
m_overRho <- as.numeric(t(a)%*%Inf_eps%*%a)
Den<- (m_overRho*meanL1^2/m2^2*var*r^2+t(a)%*%Q0_overRho+t(a)%*%P0_overRho%*%a+1)
num <-(meanL1^2/m2^2*var-2*mu^2)*r^2
@@ -402,6 +409,7 @@
Inf_eps1 <- Inf_eps*rh0
Ph <- mu^2*term%*%term1%*%invB%*%term2%*%Inf_eps1
Qh <- mu*term%*%term1%*%(-term2%*%Inf_eps1-invB%*%term2%*%Inf_eps1%*%t(B_tilde))%*%e
+ # Qh <- mu*term%*%term1%*%(-term2%*%Inf_eps1-invB%*%term2%*%Inf_eps1%*%t(B))%*%e
m <- m_overRho*rh0
if(type=="correlation"){
Modified: pkg/yuima/R/cogarchNoise.R
===================================================================
--- pkg/yuima/R/cogarchNoise.R 2015-03-01 21:23:47 UTC (rev 364)
+++ pkg/yuima/R/cogarchNoise.R 2015-03-09 09:56:52 UTC (rev 365)
@@ -65,7 +65,7 @@
auxcogarch.noise<-function(cost,b,acoeff,mu,Data,freq){
res<-StationaryMoments(cost,b,acoeff,mu)
- ExpY0<-res$ExpVar
+ ExpY0<-res$ExpStatVar
q<-length(b)
a <- e <- matrix(0,nrow=q,ncol=1)
@@ -76,13 +76,13 @@
DeltaG <- diff(Data)
squaredG <- DeltaG^2
- Process_Y <- res$ExpVar
+ Process_Y <- ExpY0
var_V<-cost + sum(acoeff*Process_Y)
delta <- 1/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}
Process_Y <- cbind(Process_Y, (expm(B*delta)%*%(Process_Y[,t-1]+e*squaredG[t-1])))
- var_V[t] <- cost + sum(acoeff*Process_Y[,t])
+ var_V[t] <- cost + sum(a*Process_Y[,t])
}
#\Delta L_{t}=\frac{\Delta G_{t}}{\sqrt{V_{t}}}.
Modified: pkg/yuima/R/simulate.R
===================================================================
--- pkg/yuima/R/simulate.R 2015-03-01 21:23:47 UTC (rev 364)
+++ pkg/yuima/R/simulate.R 2015-03-09 09:56:52 UTC (rev 365)
@@ -562,7 +562,7 @@
tavect<-t(avect)
ncolsim <- (info at q+2)
- sim <- matrix(0,n,ncolsim)
+ sim <- matrix(0,n+1,ncolsim)
par.len <- length(model at parameter@all)
if(missing(true.parameter) & par.len>0){
@@ -610,7 +610,8 @@
#Simulating jump times
- intensity <- lambda*Delta
+ #intensity <- lambda*Delta
+ intensity<-lambda
jump_time<-numeric()
jump_time[1] <- rexp(1, rate = intensity)
# In yuima this part is evaluated using function eval
@@ -618,74 +619,97 @@
Time[1] <- jump_time[1]
j <- 1
numb_jum<-numeric()
- for (i in c(1:n) ){
- numb_jum[i]<-0
- while(Time[j]<i){
- numb_jum[i]<-numb_jum[i]+1
- jump_time[j+1]<-rexp(1,rate=intensity)
- Time[j+1]<-Time[j]+jump_time[j+1]
- j<-j+1
- }
+# for (i in c(1:n) ){
+# numb_jum[i]<-0
+# while(Time[j]<i){
+# numb_jum[i]<-numb_jum[i]+1
+# jump_time[j+1]<-rexp(1,rate=intensity)
+# Time[j+1]<-Time[j]+jump_time[j+1]
+# j<-j+1
+# }
+# }
+
+ while(Time[j] < Terminal){
+ jump_time[j+1]<-rexp(1,rate=intensity)
+ Time[j+1]<-Time[j]+jump_time[j+1]
+ j<-j+1
}
- total_NumbJ <- j-1
+
+ total_NumbJ <- j
# Counting the number of jumps
- N<-matrix(1,n,1)
- N[1,1]<-numb_jum[1]
- for(i in c(2:n)){
- N[i,1]=N[i-1,1]+numb_jum[i]
- }
+# N<-matrix(1,n,1)
+# N[1,1]<-numb_jum[1]
+# for(i in c(2:n)){
+# N[i,1]=N[i-1,1]+numb_jum[i]
+# }
# Simulating the driving process
F <- suppressWarnings(parse(text=gsub("^d(.+?)\\(.+?,", "r\\1(total_NumbJ,", model at measure$df$expr, perl=TRUE)))
assign("total_NumbJ",total_NumbJ, envir=yuimaEnv)
dL<-eval(F, envir=yuimaEnv)
#dL<-rnorm(total_NumbJ,mean=0,sd=1)
- L<-matrix(1,total_NumbJ,1)
- L[1]<-dL[1]
- for(j in c(2:total_NumbJ)){
- L[j]<-L[j-1] + dL[j]
- }
+# L<-matrix(1,total_NumbJ,1)
+# L[1]<-dL[1]
+# for(j in c(2:total_NumbJ)){
+# L[j]<-L[j-1] + dL[j]
+# }
# Computing the processes V and Y at jump
V<-matrix(1,total_NumbJ,1)
Y<-matrix(1,info at q,total_NumbJ)
- Y[,1]<-matrix(1,info at q,1) #Starting point for unobservable State Space Process.
- Y[,1]<-expm(AMatrix*jump_time[1])%*%Y[,1]
- Y[,1]<-Y[,1]+(value.a0+sum(tavect*Y[,1]))*evect*dL[1]^2
+ Y[,1]<-matrix(sim[1,c(3:info at q)],info at q,1) #Starting point for unobservable State Space Process.
+
V[1,]<-value.a0+sum(tavect*Y[,1])
+ G<-matrix(1, total_NumbJ,1)
+ G[1]<-0
+
for(j in c(2:total_NumbJ)){
- Y[,j]<-Y[,j-1]
- Y[,j]<-expm(AMatrix*jump_time[j])%*%Y[,j]
- Y[,j]<-Y[,j]+(value.a0+sum(tavect*Y[,j]))*evect*dL[j]^2
+ Y[,j]<-as.numeric(expm(AMatrix*jump_time[j])%*%Y[,j-1])+(V[j-1,])*evect*dL[j-1]^2
V[j,]<-value.a0+sum(tavect*Y[,j])
- }
- # Computing the process G at jump time
- G<-matrix(1, total_NumbJ,1)
- G[1]<-sqrt(V[1])*dL[1]
- for(j in c(2:total_NumbJ)){
+ # }
+# # Computing the process G at jump time
+#
+# for(j in c(2:total_NumbJ)){
G[j]<-G[j-1]+sqrt(V[j])*dL[j]
}
- # Realizations observed at integer times
- i<-1
- while(N[i]==0){
- i <- i+1
- }
-# G_obs<-numeric()
- L_obs<-numeric()
-# V_obs<-numeric()
-# Y_obs<-matrix(0,info at q,)
- sim[c(1:(i-1)),1]<-0
- sim[c(i:n),1]<-G[N[c(i:n)]]
- L_obs[c(1:(i-1))]<-0
- L_obs[c(i:n)]<-L[N[c(i:n)]]
- for(j in c(1:(i-1))){
- sim[j,3:ncolsim]<-as.numeric(Y[,j])
- sim[j,2]<-value.a0+tavect%*%expm(AMatrix*j)%*%matrix(1,info at q,1)#Starting point for unobservable State Space Process
- }
- for(j in c(i:n)){
- sim[j,3:ncolsim]<-as.numeric(Y[,N[j]])
- sim[j,2]<-value.a0+as.numeric(tavect%*%expm(AMatrix*(j-Time[N[j]]))%*%Y[,N[j]])
- }
+
+
+ res<-approx(x=c(0,Time), y = c(0,G),
+ xout=seq(0,Terminal, by=Terminal/n),
+ method = "constant")
+ sim[,1]<-res$y
+ i<-1
+ for(j in 1:length(Time)){
+ while (i*Delta < Time[j] && i <= n){
+ sim[i+1,3:ncolsim]<-expm(AMatrix*(Time[j]-i*Delta))%*%Y[,j]
+ sim[i+1,2]<-value.a0+as.numeric(tavect%*%sim[i,3:ncolsim])
+ i<-i+1
+
+ }
+ }
+
+
+# # Realizations observed at integer times
+# i<-1
+# while(N[i]==0){
+# i <- i+1
+# }
+# # G_obs<-numeric()
+# # L_obs<-numeric()
+# # V_obs<-numeric()
+# # Y_obs<-matrix(0,info at q,)
+# sim[c(1:(i-1)),1]<-0
+# sim[c(i:n),1]<-G[N[c(i:n)]]
+# # L_obs[c(1:(i-1))]<-0
+# # L_obs[c(i:n)]<-L[N[c(i:n)]]
+# for(j in c(1:(i-1))){
+# sim[j,3:ncolsim]<-as.numeric(Y[,j])
+# sim[j,2]<-value.a0+tavect%*%expm(AMatrix*j)%*%matrix(1,info at q,1)#Starting point for unobservable State Space Process
+# }
+# for(j in c(i:n)){
+# sim[j,3:ncolsim]<-as.numeric(Y[,N[j]])
+# sim[j,2]<-value.a0+as.numeric(tavect%*%expm(AMatrix*(j-Time[N[j]]))%*%Y[,N[j]])
+# }
}
- X <- ts(sim)
+ X <- ts(sim[-1,])
Data <- setData(X,delta = Delta)
result <- setYuima(data=Data,model=yuimaCogarch at model, sampling=yuimaCogarch at sampling)
}
More information about the Yuima-commits
mailing list