[Yuima-commits] r295 - in pkg/yuima: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Mar 27 15:57:48 CET 2014
Author: lorenzo
Date: 2014-03-27 15:57:48 +0100 (Thu, 27 Mar 2014)
New Revision: 295
Removed:
pkg/yuima/R/CarmaRecovNoise.R
pkg/yuima/man/CarmaRecovNoise.Rd
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/R/qmle.R
Log:
Modified qmle.R for CARMA model
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2014-03-23 15:57:44 UTC (rev 294)
+++ pkg/yuima/DESCRIPTION 2014-03-27 14:57:48 UTC (rev 295)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project package
-Version: 1.0.5
+Version: 1.0.6
Date: 2014-03-23
Depends: methods, zoo, stats4, utils, expm
Suggests: cubature, mvtnorm
Deleted: pkg/yuima/R/CarmaRecovNoise.R
===================================================================
--- pkg/yuima/R/CarmaRecovNoise.R 2014-03-23 15:57:44 UTC (rev 294)
+++ pkg/yuima/R/CarmaRecovNoise.R 2014-03-27 14:57:48 UTC (rev 295)
@@ -1,366 +0,0 @@
-
-# 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),]
-}
-
-
-
-yuima.carma.eigen<-function(A){
- diagA<-eigen(A)
- diagA$values<-diagA$values[order(diagA$values, na.last = TRUE, decreasing = TRUE)]
- n_eigenval<-length(diagA$values)
- 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,
- matrix(diagA$values[i]^(c(1:n_eigenval)-1),n_eigenval,1))
- }
- }
- return(diagA)
-}
-
-
-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)
- e_q<-rep(0,q)
- e_q[q]<-1
- X_0<-matrix(X_0,q,1)
- e_q<-matrix(e_q,q,1)
- X<-matrix(0,q,Time)
- int<-matrix(0,q,1)
- 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){
-# # 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
- 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{
- 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))])
- }
-}
-
-
-bEvalPoly<-function(b,lambdax){
- result<-sum(b*lambdax^(c(1:length(b))-1))
- return(result)
-}
-
-aEvalPoly<-function(a,lambdax){
- p<-length(a)
- a.new<-c(1,a[c(1:(p-1))])
- pa.new<-c(p:1)*a.new
- result<-sum(pa.new*lambdax^(c(p:1)-1))
- return(result)
-}
-
-
-CarmaRecovNoise<-function(yuima, param, data=NULL,NoNeg.Noise=FALSE){
- 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)){
- observ<-yuima at data
- }else{observ<-data}
-}else{
- if(is(yuima,"yuima.carma")){
- model<-yuima
- if(is.null(data)){
- yuima.stop("Missing data")
- }
- observ<-data
- }
-}
-
- if(!is(observ,"yuima.data")){
- 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.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,
- 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)){
- 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")
- 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_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
- 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
- X_0<-rep(0,(q-1))
- }else{
- e_q<-1
- X_0<-0
- }
-
-
- 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
- #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)
-
- # 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])
- 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))
-
- 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_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-03-23 15:57:44 UTC (rev 294)
+++ pkg/yuima/R/qmle.R 2014-03-27 14:57:48 UTC (rev 295)
@@ -681,7 +681,7 @@
final_res<-new("mle", call = call, coef = coef, fullcoef = unlist(mycoef),
vcov = vcov, min = min, details = oout, minuslogl = minusquasilogl,
method = method)
- }else{
+ }else{
if( Est.Incr=="Carma.IncPar" || Est.Incr=="Carma.Inc" ){
final_res<-new("yuima.carma.qmle", call = call, coef = coef, fullcoef = unlist(mycoef),
vcov = vcov, min = min, details = oout, minuslogl = minusquasilogl,
@@ -704,7 +704,7 @@
if(!is(yuima at model,"yuima.carma")){
return(final_res)
}else {
- cat("\nStarting Estimation Increments ...\n")
+
param<-coef(final_res)
observ<-yuima at data
@@ -733,6 +733,12 @@
if (length(info at lin.par)!=0){
lin.par<-param[info at lin.par]
}
+ if(min(yuima.PhamBreton.Alg(ar.par[numb.ar:1]))>=0){
+ cat("\n Stationarity condition is satisfied...\n Starting Estimation Increments ...\n")
+ }else{
+ yuima.warn("Insert constraints in Autoregressive parameters for enforcing stationarity" )
+ cat("\n Starting Estimation Increments ...\n")
+ }
ttt<-observ at zoo.data[[1]]
tt<-index(ttt)
@@ -1512,7 +1518,7 @@
# SigMatr<-expA%*%V_inf%*%t(expA)+Qmatr
#SigMatr<-Qmatr
- #SigMatr<-V_inf
+ SigMatr<-V_inf
zc<-matrix(bvector,1,p)
loglstar <- 0
@@ -1543,6 +1549,54 @@
+yuima.PhamBreton.Alg<-function(a){
+ p<-length(a)
+ gamma<-a[p:1]
+ if(p>2){
+ gamma[p]<-a[1]
+ alpha<-matrix(NA,p,p)
+ for(j in 1:p){
+ if(is.integer(as.integer(j)/2)){
+ alpha[p,j]<-0
+ alpha[p-1,j]<-0
+ }else{
+ alpha[p,j]<-a[j]
+ alpha[p-1,j]<-a[j+1]/gamma[p]
+ }
+ }
+ for(n in (p-1):1){
+ gamma[n]<-alpha[n+1,2]-alpha[n,2]
+ for(j in 1:n-1){
+ alpha[n-1,j]<-(alpha[n+1,j+2]-alpha[n,j+2])/gamma[n]
+ }
+ alpha[n-1,n-1]<-alpha[n+1,n+1]/gamma[n]
+ }
+ gamma[1]<-alpha[2,2]
+ }
+ return(gamma)
+}
+
+#yuima.PhamBreton.Inv<-function(gamma){
+# p<-length(gamma)
+# a<-gamma[p:1]
+# if(p>2){
+# x<-polynom()
+# f0<-1*x^0
+# f1<-x
+# f2<-x*f1+gamma[1]*f0
+# for(t in 2:(p-1)){
+# f0<-f1
+# f1<-f2
+# f2<-x*f1+gamma[t]*f0
+# }
+# finpol<-f2+gamma[p]*f1
+# a <- coef(finpol)[p:1]
+# }
+# return(a)
+# }
+
+
+
#yuima.carma.loglik1<-function (y, tt, a, b, sigma)
yuima.carma.loglik1<-function (y, u, a, b, sigma,time.obs,V_inf0,p,q)
{
Deleted: pkg/yuima/man/CarmaRecovNoise.Rd
===================================================================
--- pkg/yuima/man/CarmaRecovNoise.Rd 2014-03-23 15:57:44 UTC (rev 294)
+++ pkg/yuima/man/CarmaRecovNoise.Rd 2014-03-27 14:57:48 UTC (rev 295)
@@ -1,131 +0,0 @@
-\name{CarmaRecovNoise}
-\alias{Recovering.Noise}
-\alias{Carma.Recovering}
-\alias{CarmaRecovNoise}
-\alias{Levy.Carma}
-\title{Estimation for the underlying Levy in a carma model}
-\description{Retrieve the increment of the underlying Levy for the carma(p,q) process using the approach developed in Brockwell et al.(2011)}
-\usage{
-CarmaRecovNoise(yuima, param, data=NULL, NoNeg.Noise=FALSE)
-}
-\arguments{
- \item{yuima}{a yuima object or an object of \code{\link{yuima.carma-class}}.}
- \item{param}{\code{list} of parameters for the carma.}
- \item{data}{an object of class \code{\link{yuima.data-class}} contains the observations available at uniformly spaced time. If \code{data=NULL}, the default, the 'CarmaRecovNoise' uses the data in an object of \code{\link{yuima.data-class}}.}
- \item{NoNeg.Noise}{Estimate a non-negative Levy-Driven Carma process. By default \code{NoNeg.Noise=FALSE}.}
-}
-\value{
- %\item{QL}{a real value.}
- \item{incr.Levy}{a numeric object contains the estimated increments.}
- }
-\author{The YUIMA Project Team}
-\note{
- %The function ml.ql uses the function optim internally.
- The function \code{qmle} uses the function \code{CarmaRecovNoise} for estimation of underlying Levy in the carma model.
-}
-
- \references{
- Brockwell, P., Davis, A. R. and Yang. Y. (2011)
-Estimation for Non-Negative Levy-Driven CARMA Process, \emph{Journal of Business And Economic Statistics}, \bold{29} - 2, 250-259.
-}
-
-
-\examples{
-\dontrun{
-#Ex.1: Carma(p=3, q=0) process driven by a brownian motion.
-
-mod0<-setCarma(p=3,q=0)
-
-# We fix the autoregressive and moving average parameters
-# to ensure the existence of a second order stationary solution for the process.
-
-true.parm0 <-list(a1=4,a2=4.75,a3=1.5,b0=1)
-
-# We simulate a trajectory of the Carma model.
-
-numb.sim<-1000
-samp0<-setSampling(Terminal=100,n=numb.sim)
-set.seed(100)
-incr.W<-matrix(rnorm(n=numb.sim,mean=0,sd=sqrt(100/numb.sim)),1,numb.sim)
-
-sim0<-simulate(mod0,
- true.parameter=true.parm0,
- sampling=samp0, increment.W=incr.W)
-
-#Applying the CarmaRecovNoise
-
-system.time(
- inc.Levy0<-CarmaRecovNoise(sim0,true.parm0)
-)
-
-# We compare the orginal with the estimated noise increments
-
-par(mfrow=c(1,2))
-plot(t(incr.W)[1:998],type="l", ylab="",xlab="time")
-title(main="True Brownian Motion",font.main="1")
-plot(inc.Levy0,type="l", main="Filtered Brownian Motion",font.main="1",ylab="",xlab="time")
-
-# Ex.2: carma(2,1) driven by a compound poisson
-# where jump size is normally distributed and
-# the lambda is equal to 1.
-
-mod1<-setCarma(p=2,
- q=1,
- measure=list(intensity="Lamb",df=list("dnorm(z, 0, 1)")),
- measure.type="CP")
-
-true.parm1 <-list(a1=1.39631, a2=0.05029,
- b0=1,b1=2,
- Lamb=1)
-
-# We generate a sample path.
-
-samp1<-setSampling(Terminal=100,n=200)
-set.seed(123)
-sim1<-simulate(mod1,
- true.parameter=true.parm1,
- sampling=samp1)
-
-# We estimate the parameter using qmle.
-carmaopt1 <- qmle(sim1, start=true.parm1)
-summary(carmaopt1)
-# Internally qmle uses CarmaRecovNoise. The result is in
-plot(carmaopt1)
-
-# Ex.3: Carma(p=2,q=1) with scale and location parameters
-# driven by a Compound Poisson
-# with jump size normally distributed.
-mod2<-setCarma(p=2,
- q=1,
- loc.par="mu",
- scale.par="sig",
- measure=list(intensity="Lamb",df=list("dnorm(z, 0, 1)")),
- measure.type="CP")
-
-true.parm2 <-list(a1=1.39631,
- a2=0.05029,
- b0=1,
- b1=2,
- Lamb=1,
- mu=0.5,
- sig=0.23)
-# We simulate the sample path
-set.seed(123)
-sim2<-simulate(mod2,
- true.parameter=true.parm2,
- sampling=samp1)
-
-# We estimate the Carma and we plot the underlying noise.
-
-carmaopt2 <- qmle(sim2, start=true.parm2)
-summary(carmaopt2)
-
-# Increments estimated by CarmaRecovNoise
-plot(carmaopt2)
-}
-}
-
-
-% Add one or more standard keywords, see file 'KEYWORDS' in the
-% R documentation directory.
-\keyword{ts}
More information about the Yuima-commits
mailing list