[Yuima-commits] r268 - pkg/yuima/R

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Tue Jan 21 22:49:23 CET 2014


Author: lorenzo
Date: 2014-01-21 22:49:22 +0100 (Tue, 21 Jan 2014)
New Revision: 268

Modified:
   pkg/yuima/R/CarmaRecovNoise.R
   pkg/yuima/R/qmle.R
Log:
Bugs Removed

Modified: pkg/yuima/R/CarmaRecovNoise.R
===================================================================
--- pkg/yuima/R/CarmaRecovNoise.R	2014-01-13 21:07:48 UTC (rev 267)
+++ pkg/yuima/R/CarmaRecovNoise.R	2014-01-21 21:49:22 UTC (rev 268)
@@ -4,6 +4,21 @@
 # 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),]
+}
+
+
+
 yuima.carma.eigen<-function(A){
   diagA<-eigen(A)
   diagA$values<-diagA$values[order(diagA$values, na.last = TRUE, decreasing = TRUE)]
@@ -19,7 +34,7 @@
 }
 
 
-StateVarX<-function(y,tt,X_0,B){
+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) 
   q<-length(X_0)
@@ -29,39 +44,92 @@
   e_q<-matrix(e_q,q,1)
   X<-matrix(0,q,Time)
   int<-matrix(0,q,1)
-  for (i in c(2:Time)){
-    int<-int+(expm(B*(tt[i]-tt[(i-1)]))%*%(e_q*y[i]))*(tt[i]-tt[(i-1)])
-    X[i]<-as.matrix(expm(B*tt[i])%*%X_0+int)
+  for (i in c(3:Time)){
+    if(discr.eul==FALSE){
+      int<-int+(expm(B*(tt[i]-tt[(i-1)]))%*%(e_q*y[i-1]))*(tt[i]-tt[(i-1)])
+      X[,i]<-as.matrix(expm(B*tt[i])%*%X_0+int)
+    }else{
+      X[,i]<-X[,i-1]+(B%*%X[,i-1])*(tt[i]-tt[(i-1)])+(e_q*y[i-1])*(tt[i]-tt[(i-1)])
+    }
   }
   return(X)
 }
 
 
-StateVarXp<-function(y,X_q,tt,B,q,p){
+# StateVarXp<-function(y,X_q,tt,B,q,p){
+#   # 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)
+#   }else{
+#     MatrD<-as.matrix(diagMatB$values)
+#   }
+#   MatrR<-diagMatB$vectors
+#   idx.r<-c(q:(p-1))
+#   elem.X <-length(idx.r)
+#   YMatr<-matrix(0,q,length(y))
+#   YMatr[q,]<-y
+#   OutherX<-matrix(NA,elem.X,length(y))
+#   # 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,]  
+#   }
+#   X.StatVar<-rbind(X_q,OutherX)
+#   return(X.StatVar)
+# }
+
+
+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 
   # see Brockwell, Davis and Yang 2011
-  
-  diagMatB<-yuima.carma.eigen(B)
-  if(length(diagMatB$values)>1){
-    MatrD<-diag(diagMatB$values)
+  if(nume.Der==FALSE){
+    yuima.warn("We need to develop this part in the future for gaining speed")
+    return(NULL)
+#     diagMatB<-yuima.carma.eigen(B)
+#     if(length(diagMatB$values)>1){
+#       MatrD<-diag(diagMatB$values)
+#     }else{
+#       MatrD<-as.matrix(diagMatB$values)
+#     }
+#     MatrR<-diagMatB$vectors
+#     idx.r<-c(q:(p-1))
+#     elem.X <-length(idx.r)
+#     YMatr<-matrix(0,q,length(y))
+#     YMatr[q,]<-y
+#     OutherX<-matrix(NA,elem.X,length(y))
+#     # 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,]  
+#     }
+#     X.StatVar<-rbind(X_q,OutherX)
+#     return(X.StatVar)
   }else{
-    MatrD<-as.matrix(diagMatB$values)
+    X_q0<-as.numeric(X_q[1,])
+    qlen<-length(X_q0)
+    X.StatVar<-matrix(0,p,qlen)
+    X.StatVar[1,]<-X_q0[1:length(X_q0)]
+    diffStatVar<-diff(X_q0)
+    difftime<-diff(tt)[1]
+    for(i in 2:p){
+      in.count<-(p-i+1)
+      fin.count<-length(diffStatVar)
+      dummyDer<-(diffStatVar/difftime)
+      X.StatVar[i,c(1:length(dummyDer))]<-(dummyDer)
+      diffStatVar<-diff(dummyDer)
+#       if(in.count-1>0){
+#         difftime<-difftime[c((in.count-1):fin.count)]
+#       }else{difftime<-difftime[c((in.count):fin.count)]}
+    }
+    X.StatVar[c(1:dim(X_q)[1]),]<-X_q
+    return(X.StatVar[,c(1:length(dummyDer))])
   }
-  MatrR<-diagMatB$vectors
-  idx.r<-c(q:(p-1))
-  elem.X <-length(idx.r)
-  YMatr<-matrix(0,q,length(y))
-  YMatr[q,]<-y
-  OutherX<-matrix(NA,q,length(y))
-  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,]  
-  }
-  X.StatVar<-rbind(X_q,OutherX)
-  return(X.StatVar)
 }
 
+
 bEvalPoly<-function(b,lambdax){
   result<-sum(b*lambdax^(c(1:length(b))-1))
   return(result)
@@ -76,7 +144,7 @@
 }
 
 
-CarmaRecovNoise<-function(yuima, param, data=NULL){
+CarmaRecovNoise<-function(yuima, param, data=NULL,NoNeg.Noise=FALSE){
   if( missing(param) ) 
     yuima.stop("Parameter values are missing.")
   
@@ -137,21 +205,39 @@
   tt<-index(ttt)
   y<-coredata(ttt)
   
-  levy<-yuima.CarmaRecovNoise(y,tt,ar.par,ma.par, loc.par, scale.par, lin.par)
-  inc.levy<-diff(t(levy))
+  levy<-yuima.CarmaRecovNoise(y,tt,ar.par,ma.par, loc.par, scale.par, lin.par,NoNeg.Noise)
+  inc.levy<-diff(as.numeric(levy))
   return(inc.levy)
 }
 
 
 
-yuima.CarmaRecovNoise<-function(y,tt,ar.par,ma.par, loc.par=NULL, scale.par=NULL, lin.par=NULL){
-    if(!is.null(loc.par)){
-      yuima.warn("the loc.par will be implemented as soon as possible")
-      return(NULL)
+yuima.CarmaRecovNoise<-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")
+        #mean.y<-tail(ma.par,n=1)/ar.par[1]*scale.par
+        mean.y<-mean(y)
+        #y<-y-mean.y
+        mean.L1<-mean.y/tail(ma.par,n=1)*ar.par[1]/scale.par
+      }
+    }  
     if(!is.null(scale.par)){
-      yuima.warn("the scale.par will be implemented as soon as possible")
-      return(NULL)
+      ma.parAux<-ma.par*scale.par
+      names(ma.parAux)<-names(ma.par)
+      ma.par<-ma.parAux
+#       yuima.warn("the scale.par will be implemented as soon as possible")
+#       return(NULL)
     }
     if(!is.null(lin.par)){
       yuima.warn("the lin.par will be implemented as soon as possible")
@@ -161,17 +247,34 @@
     p<-length(ar.par)
     q<-length(ma.par)
     
-    if(q==1){
-      yuima.warn("the car(p) process will be implemented as soon as possible")
-      return(NULL)
-    }else{
     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_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)){          
+          dummyDerY<-diffY/diffX
+          X.StVa[i,c(1:length(dummyDerY))]<-(dummyDerY)
+          diffY<-diff(dummyDerY)
+        
+        }
+        X.StVa<-X.StVa[,c(1:length(dummyDerY))]
+      }
+      #yuima.warn("the car(p) process will be implemented as soon as possible")
+      #return(NULL)
+    }else{
     newma.par<-ma.par/b_q
     # We build the matrix B which is necessary for building eq 5.2
-    ynew<-y/b_q
     B<-MatrixA(newma.par[c(1:(q-1))])
     diagB<-yuima.carma.eigen(B)
     
@@ -185,20 +288,23 @@
     }
     
     
- 
-    X_q<-StateVarX(ynew,tt,X_0,B)
+    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)
-    
+    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])
@@ -209,23 +315,52 @@
     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
-    # yver<-Y_CVS[1,]+Y_CVS[2,]
+    # yver1<-Y_CVS[1,]+Y_CVS[2,]+Y_CVS[3,]
     
-    # plot(yver)
+    # plot(yver1)
     
     # plot(y)
+    # Prop 2 Verified even in the case of q=0
     
-    
     idx.r<-match(0,Im(diagA$values))
-    lambda.r<-diagA$values[idx.r]
+    lambda.r<-Re(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 
+#     }
+    if(nume.Der==TRUE){
+      tt<-tt[p:length(tt)]
+      # CHECK HERE 15/01 
+    }
+    
     lev.und<-matrix(0,1,length(tt))
-    derA<-aEvalPoly(ar.par[c(p:1)],lambda.r)
+      
+    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)]
+                #                                                            +(matrix(c(0,mean.L1),2,(dim(X.StVa)[2]-1)))/diff(tt)[1]
+                                                                            )*diff(tt)[1]
+      #+(matrix(c(0,mean.L1),2,(dim(X.StVa)[2]-1)))
+      lev.und<-as.matrix(cumsum(as.numeric(incr.vect[dim(X.StVa)[1],])),1,length(tt))
+    }else{
       for(t in c(2:length(tt))){
-        int<-int+y[t]*(tt[t]-tt[t-1])
+        int<-int+Y_CVS[idx.r,t-1]*(tt[t]-tt[t-1])
         lev.und[,t]<-derA/BinLambda[idx.r]*(Y_CVS[idx.r,t]-Y_CVS[idx.r,1]-lambda.r*int)
       }
     }
+#     if(NoNeg.Noise==TRUE){
+#       if(length(ar.par)==length(ma.par)){
+#         if(Alternative==TRUE){
+#           mean.Ltt<-mean.L1*tt[c(2:length(tt))]
+#         }else{mean.Ltt<-mean.L1*tt}
+#         lev.und<-lev.und+mean.Ltt
+#       
+#       }
+#     }
   return(lev.und)
 }
 
+

Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R	2014-01-13 21:07:48 UTC (rev 267)
+++ pkg/yuima/R/qmle.R	2014-01-21 21:49:22 UTC (rev 268)
@@ -75,7 +75,9 @@
 
 qmle <- function(yuima, start, method="BFGS", fixed = list(), print=FALSE, 
  lower, upper, joint=FALSE, ...){
-
+  if(is(yuima at model, "yuima.carma")){
+    NoNeg.Noise<-FALSE
+  }
   if(is(yuima at model, "yuima.carma")&& length(yuima at model@info at scale.par)!=0){
     method<-"L-BFGS-B"
   }
@@ -118,15 +120,24 @@
       tmp <- regexpr("\\(", yuima at model@measure$df$exp)[1]
       measurefunc <- substring(yuima at model@measure$df$exp, 1, tmp-1)
       if(!is.na(match(measurefunc,CPlist))){
-        yuima.warn("carma(p,q): the qmle for a carma(p,q) driven by a Compound Poisson with no-negative random size will be implemented as soon as possible")
-        return(NULL)
+        yuima.warn("carma(p,q): the qmle for a carma(p,q) driven by a Compound Poisson with no-negative random size")
+        NoNeg.Noise<-TRUE 
+        # we need to add a new parameter for the mean of the Noise
+        if((yuima at model@info at q+1)==(yuima at model@info at q+1)){
+          start[["mean.noise"]]<-1
+        }
+  #      return(NULL)
       }
       if(yuima at model@measure.type=="code"){
         tmp <- regexpr("\\(", yuima at model@measure$df$exp)[1]
         measurefunc <- substring(yuima at model@measure$df$exp, 1, tmp-1)
         if(!is.na(match(measurefunc,codelist))){
-          yuima.warn("carma(p,q): the qmle for a carma(p,q) driven by a Compound Poisson with no-negative random size will be implemented as soon as possible")
-          return(NULL)
+          yuima.warn("carma(p,q): the qmle for a carma(p,q) driven by a non-Negative Levy  will be implemented as soon as possible")
+          NoNeg.Noise<-TRUE
+          if((yuima at model@info at q+1)==(yuima at model@info at q+1)){
+            start[[mean.noise]]<-1
+          }
+          #return(NULL)
         }      
       }
     }
@@ -204,11 +215,32 @@
     fullcoef<-c(fullcoef, yuima at model@info at loc.par)
     
   }
+
+  if(is(yuima at model,"yuima.carma") && (NoNeg.Noise==TRUE)){
+    if((yuima at model@info at q+1)==yuima at model@info at p){
+      mean.noise<-"mean.noise"
+      fullcoef<-c(fullcoef, mean.noise)
+    }
+  }
   
+  if(is(yuima at model,"yuima.carma") && (length(measure.par)>0)){
+   fullcoef<-c(fullcoef, measure.par)
+  }
+
 	npar <- length(fullcoef)
 	
-	fixed.par <- names(fixed)
-
+  
+	fixed.par <- names(fixed) # We use Fixed.par when we consider a Carma with scale parameter
+  if(is(yuima at model,"yuima.carma") && (length(measure.par)>0)){
+    
+    for( j in c(1:length(measure.par))){
+      if(is.na(match(measure.par[j],names(fixed)))){
+        fixed.par <- c(fixed.par,measure.par[j])
+        fixed[measure.par[j]]<-start[measure.par[j]]
+      }
+    }
+    
+  }
 	if (any(!(fixed.par %in% fullcoef))) 
 	 yuima.stop("Some named arguments in 'fixed' are not arguments to the supplied yuima model")
 	
@@ -228,7 +260,7 @@
 	#01/01
   if(is(yuima at model, "yuima.carma")){
    # if(length(yuima at model@info at scale.par)!=0){
-      idx.xinit <- as.integer(na.omit(match(xinit.par,nm)))
+      idx.xinit <- as.integer(na.omit(match(xinit.par,nm)))# We need to add idx if NoNeg.Noise is TRUE
     #}
   }
   
@@ -254,6 +286,7 @@
 	}
 	upper <- tmpupper
 
+  
 	 
 	 
 	 
@@ -307,7 +340,9 @@
 	 HaveDiffHess <- FALSE
 	 if(length(start)){
 		if(JointOptim){ ### joint optimization
-			if(length(start)>1){ #Â?multidimensional optim
+			if(length(start)>1){ #Â?multidimensional optim # Adjust lower for no negative Noise
+        if(is(yuima at model,"yuima.carma") && (NoNeg.Noise==TRUE))
+          if(mean.noise %in% names(lower)){lower[mean.noise]<-10^-7}
 				oout <- optim(start, fj, method = method, hessian = TRUE, lower=lower, upper=upper)
 				HESS <- oout$hessian
 				if(is(yuima at model,"yuima.carma") && length(yuima at model@info at scale.par)!=0){
@@ -316,6 +351,18 @@
 				  HESS<-HESS[-idx.b0,]
 				  HESS<-HESS[,-idx.b0]
 				}
+				if(is(yuima at model,"yuima.carma") && length(yuima at model@parameter at measure)!=0){
+				  for(i in c(1:length(fixed.par))){
+            indx.fixed<-match(fixed.par[i],rownames(HESS))
+            HESS<-HESS[-indx.fixed,]
+				    HESS<-HESS[,-indx.fixed]
+				  }
+          if(is(yuima at model,"yuima.carma") && (NoNeg.Noise==TRUE)){
+            idx.noise<-(match(mean.noise,rownames(HESS)))
+            HESS<-HESS[-idx.noise,]
+            HESS<-HESS[,-idx.noise]
+          }
+				}
 				HaveDriftHess <- TRUE
 				HaveDiffHess <- TRUE
 			} else { ### one dimensional optim
@@ -394,7 +441,13 @@
 			mydots$lower <- unlist( lower[ nm[idx.drift] ])
 			if(length(mydots$par)>1){
 			  if(is(yuima at model, "yuima.carma")){
-			    if(length(yuima at model@info at scale.par)!=0){
+			    if(NoNeg.Noise==TRUE){
+            if((yuima at model@info at q+1)==yuima at model@info at p){
+              mydots$lower[names(start["NoNeg.Noise"])]<-10^(-7)
+            }
+             
+			    }
+          if(length(yuima at model@info at scale.par)!=0){
             name_b0<-paste0(yuima at model@info at ma.par,"0",collapse="")
             index_b0<-match(name_b0,nm)
   			    mydots$lower[index_b0]<-1
@@ -462,10 +515,15 @@
 	 coef <- oout$par
 	 control=list()
 	 par <- coef
-  if(!is(yuima at model,"yuima.carma")){
-	  names(par) <- c(diff.par, drift.par)
-       nm <- c(diff.par, drift.par)
+  #if(!is(yuima at model,"yuima.carma")){
+	  names(par) <- unique(c(diff.par, drift.par))
+       nm <- unique(c(diff.par, drift.par))
+  if(is(yuima at model,"yuima.carma") && (NoNeg.Noise==TRUE)){
+      nm <-c(nm,mean.noise)
+      nm <-c(nm,measure.par)
+      nm<-unique(nm)
   }
+  #}
 #	 print(par)
 #	 print(coef)
 	 conDrift <- list(trace = 5, fnscale = 1, 
@@ -601,8 +659,9 @@
     ttt<-observ at zoo.data[[1]]
     tt<-index(ttt)
     y<-coredata(ttt)
+    if(NoNeg.Noise==TRUE && (info at p==(info at q+1))){final_res at coef[mean.noise]<-mean(y)/tail(ma.par,n=1)*ar.par[1]}
     
-    levy<-yuima.CarmaRecovNoise(y,tt,ar.par,ma.par, loc.par, scale.par, lin.par)
+    levy<-yuima.CarmaRecovNoise(y,tt,ar.par,ma.par, loc.par, scale.par, lin.par, NoNeg.Noise)
     inc.levy<-NULL
     if (!is.null(levy)){
       inc.levy<-diff(t(levy))
@@ -642,8 +701,14 @@
 	if(is(yuima at model, "yuima.carma") && length(yuima at model@info at lin.par)==0
 	   && length(yuima at model@parameter at jump)!=0){
 	  diff.par<-yuima at model@parameter at jump
+	  measure.par<-yuima at model@parameter at measure
 	}
 	
+	if(is(yuima at model, "yuima.carma") && length(yuima at model@info at lin.par)==0
+	   && length(yuima at model@parameter at measure)!=0){
+	  measure.par<-yuima at model@parameter at measure
+	}
+  
 	# 24/12
 	if(is(yuima at model, "yuima.carma") && length(yuima at model@info at lin.par)>0  ){
 	  yuima.warn("carma(p,q): the case of lin.par will be implemented as soon as")
@@ -683,6 +748,17 @@
 	        fullcoef<-c(fullcoef, xinit.par)
 				}
   
+  	if(is(yuima at model,"yuima.carma") && (length(yuima at model@parameter at measure)!=0)){
+  	  fullcoef<-c(fullcoef, measure.par)
+  	}
+  if(is(yuima at model,"yuima.carma")){
+  	if("mean.noise" %in% names(param)){
+      mean.noise<-"mean.noise"
+  	  fullcoef<-c(fullcoef, mean.noise)
+      NoNeg.Noise<-TRUE
+  	}
+  }
+  
 #	fullcoef <- c(diff.par, drift.par)
 	npar <- length(fullcoef)
 #	cat("\nparam\n")
@@ -773,6 +849,26 @@
           } 
        sigma<-tail(param,1)
          }else {sigma<-1}
+         NoNeg.Noise<-FALSE
+         if(is(yuima at model,"yuima.carma")){
+           if("mean.noise" %in% names(param)){
+             
+             NoNeg.Noise<-TRUE
+           }
+         }
+         if(NoNeg.Noise==TRUE){
+           if (length(b)==p){
+             #mean.noise<-param[mean.noise]
+             # Be useful for carma driven by a no negative levy process
+             mean.y<-mean(y)
+             #mean.y<-mean.noise*tail(b,n=1)/tail(a,n=1)*sigma
+             #param[mean.noise]<-mean.y/(tail(b,n=1)/tail(a,n=1)*sigma)
+           }else{
+             mean.y<-0
+           }
+           y<-y-mean.y
+         }
+         
        strLog<-yuima.carma.loglik1(y, tt, a, b, sigma)
        }else{
          # 01/01
@@ -795,13 +891,19 @@
          } else{sigma <- 1}
          loc.par <- yuima at model@info at loc.par
          mu <- param[loc.par]
+# Lines 883:840 work if we have a no negative noise
+        if(is(yuima at model,"yuima.carma")&&(NoNeg.Noise==TRUE)){
+           if (length(b)==p){
+             mean.noise<-param[mean.noise]
+           # Be useful for carma driven by levy process
+             mean.y<-mean.noise*tail(b,n=1)/tail(a,n=1)*sigma
+           }else{
+             mean.y<-0
+           }
+           y<-y-mean.y
+        }
          
-#          if (length(b)==p){
-#          # Be useful for carma driven by levy process
-#            mean.y<-mu*tail(b,n=1)/tail(a,n=1)*sigma
-#          }else{
-#            mean.y<-0
-#          }
+         
          y.start <- y-mu
          
          strLog<-yuima.carma.loglik1(y.start, tt, a, b, sigma)



More information about the Yuima-commits mailing list