[Yuima-commits] r615 - pkg/yuima/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jun 6 16:22:40 CEST 2017
Author: lorenzo
Date: 2017-06-06 16:22:40 +0200 (Tue, 06 Jun 2017)
New Revision: 615
Modified:
pkg/yuima/R/CarmaNoise.R
Log:
Fixed Bugs CarmaNoise.R
Modified: pkg/yuima/R/CarmaNoise.R
===================================================================
--- pkg/yuima/R/CarmaNoise.R 2017-05-28 12:22:34 UTC (rev 614)
+++ pkg/yuima/R/CarmaNoise.R 2017-06-06 14:22:40 UTC (rev 615)
@@ -1,19 +1,19 @@
-# In this file we develop the procedure described in Brockwell, Davis and Yang (2012)
-# for estimation of the underlying noise once the parameters of carma(p,q) are previously
+# In this file we develop the procedure described in Brockwell, Davis and Yang (2012)
+# for estimation of the underlying noise once the parameters of carma(p,q) are previously
# obtained
-
+
yuima.eigen2arparam<-function(lambda){
MatrixLamb<-diag(lambda)
-
+
matrixR<-matrix(NA,length(lambda),length(lambda))
for(i in c(1:length(lambda))){
matrixR[,i]<-lambda[i]^c(0:(length(lambda)-1))
}
-
+
AMatrix<-matrixR%*%MatrixLamb%*%solve(matrixR)
-
+
acoeff<-AMatrix[length(lambda),]
}
@@ -26,7 +26,7 @@
diagA$vectors<-matrix(diagA$values[1]^(c(1:n_eigenval)-1),n_eigenval,1)
if(n_eigenval>=2){
for (i in 2:n_eigenval){
- diagA$vectors<-cbind(diagA$vectors,
+ diagA$vectors<-cbind(diagA$vectors,
matrix(diagA$values[i]^(c(1:n_eigenval)-1),n_eigenval,1))
}
}
@@ -36,7 +36,7 @@
StateVarX<-function(y,tt,X_0,B,discr.eul){
# The code obtains the first q-1 state variable using eq 5.1 in Brockwell, Davis and Yang 2011
- Time<-length(tt)
+ Time<-length(tt)
q<-length(X_0)
e_q<-rep(0,q)
e_q[q]<-1
@@ -57,9 +57,9 @@
# StateVarXp<-function(y,X_q,tt,B,q,p){
-# # The code computes the state variable X using the derivatives of eq 5.2
+# # The code computes the state variable X using the derivatives of eq 5.2
# # see Brockwell, Davis and Yang 2011
-#
+#
# diagMatB<-yuima.carma.eigen(B)
# if(length(diagMatB$values)>1){
# MatrD<-diag(diagMatB$values)
@@ -75,7 +75,7 @@
# # OutherXalt<-matrix(NA,q,length(y))
# for(i in 1:elem.X){
# OutherX[i,]<-((MatrR%*%MatrD^(idx.r[i])%*%solve(MatrR))%*%X_q+
-# (MatrR%*%MatrD^(idx.r[i]-1)%*%solve(MatrR))%*%YMatr)[1,]
+# (MatrR%*%MatrD^(idx.r[i]-1)%*%solve(MatrR))%*%YMatr)[1,]
# }
# X.StatVar<-rbind(X_q,OutherX)
# return(X.StatVar)
@@ -83,7 +83,7 @@
StateVarXp<-function(y,X_q,tt,B,q,p,nume.Der){
- # The code computes the state variable X using the derivatives of eq 5.2
+ # The code computes the state variable X using the derivatives of eq 5.2
# see Brockwell, Davis and Yang 2011
if(nume.Der==FALSE){
yuima.warn("We need to develop this part in the future for gaining speed")
@@ -103,7 +103,7 @@
# # OutherXalt<-matrix(NA,q,length(y))
# for(i in 1:elem.X){
# OutherX[i,]<-((MatrR%*%MatrD^(idx.r[i])%*%solve(MatrR))%*%X_q+
-# (MatrR%*%MatrD^(idx.r[i]-1)%*%solve(MatrR))%*%YMatr)[1,]
+# (MatrR%*%MatrD^(idx.r[i]-1)%*%solve(MatrR))%*%YMatr)[1,]
# }
# X.StatVar<-rbind(X_q,OutherX)
# return(X.StatVar)
@@ -145,16 +145,16 @@
CarmaNoise<-function(yuima, param, data=NULL,NoNeg.Noise=FALSE){
- if( missing(param) )
+ if( missing(param) )
yuima.stop("Parameter values are missing.")
-
+
if(!is.list(param))
yuima.stop("Argument 'param' must be of list type.")
-
+
vect.param<-as.numeric(param)
name.param<-names(param)
names(vect.param)<-name.param
-
+
if(is(yuima,"yuima")){
model<-yuima at model
if(is.null(data)){
@@ -171,58 +171,58 @@
}
if(!is(observ,"yuima.data")){
- yuima.stop("Data must be an object of class yuima.data-class")
+ yuima.stop("Data must be an object of class yuima.data-class")
}
-
+
info<-model at info
-
+
numb.ar<-info at p
name.ar<-paste(info at ar.par,c(numb.ar:1),sep="")
ar.par<-vect.param[name.ar]
-
+
numb.ma<-info at q
name.ma<-paste(info at ma.par,c(0:numb.ma),sep="")
ma.par<-vect.param[name.ma]
-
+
loc.par=NULL
if (length(info at loc.par)!=0){
loc.par<-vect.param[info at loc.par]
}
-
+
scale.par=NULL
if (length(info at scale.par)!=0){
scale.par<-vect.param[info at scale.par]
}
-
+
lin.par=NULL
if (length(info at lin.par)!=0){
lin.par<-vect.param[info at lin.par]
}
-
-
-
+
+
+
ttt<-observ at zoo.data[[1]]
tt<-index(ttt)
y<-coredata(ttt)
-
+
levy<-yuima.CarmaNoise(y,tt,ar.par,ma.par, loc.par, scale.par, lin.par,NoNeg.Noise)
inc.levy<-diff(as.numeric(levy))
- return(inc.levy[-1]) #We start to compute the increments from the second observation.
+ return(inc.levy[-1]) #We start to compute the increments from the second observation.
}
-yuima.CarmaNoise<-function(y,tt,ar.par,ma.par,
- loc.par=NULL,
- scale.par=NULL,
- lin.par=NULL,
+yuima.CarmaNoise<-function(y,tt,ar.par,ma.par,
+ loc.par=NULL,
+ scale.par=NULL,
+ lin.par=NULL,
NoNeg.Noise=FALSE){
-
+
if(!is.null(loc.par)){
y<-y-loc.par
# yuima.warn("the loc.par will be implemented as soon as possible")
# return(NULL)
- }
+ }
if(NoNeg.Noise==TRUE){
if(length(ar.par)==length(ma.par)){
yuima.warn("The case with no negative jump needs to test deeply")
@@ -231,7 +231,7 @@
#y<-y-mean.y
mean.L1<-mean.y/tail(ma.par,n=1)*ar.par[1]/scale.par
}
- }
+ }
if(!is.null(scale.par)){
ma.parAux<-ma.par*scale.par
names(ma.parAux)<-names(ma.par)
@@ -243,30 +243,30 @@
yuima.warn("the lin.par will be implemented as soon as possible")
return(NULL)
}
-
+
p<-length(ar.par)
q<-length(ma.par)
-
+
A<-MatrixA(ar.par[c(p:1)])
b_q<-tail(ma.par,n=1)
ynew<-y/b_q
-
+
if(q==1){
yuima.warn("The Derivatives for recovering noise have been performed numerically.")
nume.Der<-TRUE
X_qtot<-ynew
- X.StVa<-matrix(0,p,length(ynew))
+ X.StVa<-matrix(0,p,length(ynew))
X_q<-X_qtot[c(1:length(X_qtot))]
-
+
X.StVa[1,]<-X_q
if(p>1){
diffY<-diff(ynew)
diffX<-diff(tt)[1]
- for (i in c(2:p)){
+ for (i in c(2:p)){
dummyDerY<-diffY/diffX
X.StVa[i,c(1:length(dummyDerY))]<-(dummyDerY)
diffY<-diff(dummyDerY)
-
+
}
X.StVa<-X.StVa[,c(1:length(dummyDerY))]
}
@@ -277,7 +277,7 @@
# We build the matrix B which is necessary for building eq 5.2
B<-MatrixA(newma.par[c(1:(q-1))])
diagB<-yuima.carma.eigen(B)
-
+
e_q<-rep(0,(q-1))
if((q-1)>0){
e_q[(q-1)]<-1
@@ -286,42 +286,47 @@
e_q<-1
X_0<-0
}
-
-
- discr.eul<-TRUE
+
+
+ discr.eul<-TRUE
# We use the Euler discretization of eq 5.1 in Brockwell, Davis and Yang
X_q<-StateVarX(ynew,tt,X_0,B,discr.eul)
-
-
+
+
nume.Der<-TRUE
#plot(t(X_q))
- X.StVa<-StateVarXp(ynew,X_q,tt,B,q-1,p,nume.Der) #Checked once the numerical derivatives have been used
+ X.StVa<-StateVarXp(ynew,X_q,tt,B,q-1,p,nume.Der) #Checked once the numerical derivatives have been used
#plot(y)
}
-
+
diagA<-yuima.carma.eigen(A)
-
-
-
-
+
+
+
+
BinLambda<-rep(NA,length(ma.par))
for(i in c(1:length(diagA$values))){
BinLambda[i]<-bEvalPoly(ma.par,diagA$values[i])
}
- MatrBLam<-diag(BinLambda)
-
+ if(length(BinLambda)==1){
+ MatrBLam<-as.matrix(BinLambda)
+ }else{
+ MatrBLam<-diag(BinLambda)
+ }
+
+
# We get the Canonical Vector Space in eq 2.17
Y_CVS<-MatrBLam%*%solve(diagA$vectors)%*%X.StVa #Canonical Vector Space
-
+
# We verify the prop 2 in the paper "Estimation for Non-Negative Levy Driven CARMA process
# yver1<-Y_CVS[1,]+Y_CVS[2,]+Y_CVS[3,]
-
+
# plot(yver1)
-
+
# plot(y)
# Prop 2 Verified even in the case of q=0
-
+
idx.r<-match(0,Im(diagA$values))
lambda.r<-Re(diagA$values[idx.r])
@@ -331,19 +336,19 @@
lambda.r<-diagA$values[idx.r]
}
int<-0
-
+
derA<-aEvalPoly(ar.par[c(p:1)],lambda.r)
# if(q==1){
# tt<-tt[p:length(tt)]
-# # CHECK HERE 15/01
+# # CHECK HERE 15/01
# }
if(nume.Der==TRUE){
tt<-tt[p:length(tt)]
- # CHECK HERE 15/01
+ # CHECK HERE 15/01
}
-
+
lev.und<-matrix(0,1,length(tt))
-
+
Alternative<-FALSE
if(Alternative==TRUE){
incr.vect<-(X.StVa[,2:dim(X.StVa)[2]]-X.StVa[,1:(dim(X.StVa)[2]-1)])-(A%*%X.StVa[,1:(dim(X.StVa)[2]-1)]
@@ -363,7 +368,7 @@
# mean.Ltt<-mean.L1*tt[c(2:length(tt))]
# }else{mean.Ltt<-mean.L1*tt}
# lev.und<-lev.und+mean.Ltt
-#
+#
# }
# }
return(Re(lev.und))
More information about the Yuima-commits
mailing list