[Yuima-commits] r789 - in pkg/yuima: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Feb 26 18:18:55 CET 2022
Author: lorenzo
Date: 2022-02-26 18:18:54 +0100 (Sat, 26 Feb 2022)
New Revision: 789
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/R/qmleLevy.R
Log:
Added New qmleLevy.R
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2022-01-28 03:51:32 UTC (rev 788)
+++ pkg/yuima/DESCRIPTION 2022-02-26 17:18:54 UTC (rev 789)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project Package for SDEs
-Version: 1.15.3
+Version: 1.15.4
Depends: R(>= 2.10.0), methods, zoo, stats4, utils, expm, cubature,
mvtnorm
Imports: Rcpp (>= 0.12.1), boot (>= 1.3-2), glassoFast,
Modified: pkg/yuima/R/qmleLevy.R
===================================================================
--- pkg/yuima/R/qmleLevy.R 2022-01-28 03:51:32 UTC (rev 788)
+++ pkg/yuima/R/qmleLevy.R 2022-02-26 17:18:54 UTC (rev 789)
@@ -299,6 +299,11 @@
myGamhat[1:dim(DiffHessian)[1],1:dim(DiffHessian)[2]]<-DiffHessian
myGamhat[dim(DiffHessian)[1]+1:dim(DriftHessian)[1],dim(DiffHessian)[1]+1:dim(DriftHessian)[2]]<-DriftHessian#293
myGamhat<-myGamhat/oldyuima at sampling@Terminal
+ }else{
+ myGamhat <- matrix(0,length(coef),length(coef))
+ myGamhat[1:dim(DiffHessian)[1],1:dim(DiffHessian)[2]]<-DiffHessian*oldyuima at sampling@Terminal/oldyuima at sampling@n
+ myGamhat[dim(DiffHessian)[1]+1:dim(DriftHessian)[1],dim(DiffHessian)[1]+1:dim(DriftHessian)[2]]<-DriftHessian#293
+ myGamhat<-myGamhat/oldyuima at sampling@Terminal
}
}else{
yuima.stop("third estimation may be theoretical invalid under the presence of an overlapping parameter.")
@@ -688,20 +693,350 @@
MatSigmaHat1<- rbind(cbind(MatSigmaHat,t(dum)),cbind(dum,SigmaEta))
InvIn<- solve(I_n)
res at vcov <-InvIn%*%MatSigmaHat1%*%t(InvIn)/Ter
+ colnames(res at vcov)<-names(res at fullcoef)
+ #res at vcov<-rbind(res at vcov,matrix(NA,nrow=length(esti$par),ncol=dim(res at vcov)[2]))
+ rownames(res at vcov)<-names(res at fullcoef)
+ res at min<-c(res at min,esti$value)
+ res at nobs<-c(res at nobs,length(Incr.Lev at zoo.data[[1]]))
+ result<- new("yuima.qmleLevy.incr",Incr.Lev=Incr.Lev,
+ minusloglLevy = minusloglik,logL.Incr=-esti$value,
+ Data = yuima at data, yuima=res, Levydetails=esti)
}else{
+ Tn<-yuima at sampling@Terminal
+ HessianEta_divTn <- HessianEta/Tn
+ if(!aggregation){
+ yuima.stop("da fare")
+ }
+
res at vcov<-cbind(res at vcov,matrix(NA,ncol=length(esti$par),nrow=dim(res at vcov)[1]))
res at vcov<-rbind(res at vcov,matrix(NA,nrow=length(esti$par),ncol=dim(res at vcov)[2]))
+ colnames(res at vcov)<-names(res at fullcoef)
+ #res at vcov<-rbind(res at vcov,matrix(NA,nrow=length(esti$par),ncol=dim(res at vcov)[2]))
+ rownames(res at vcov)<-names(res at fullcoef)
+ res at min<-c(res at min,esti$value)
+ res at nobs<-c(res at nobs,length(Incr.Lev at zoo.data[[1]]))
+
+ result<- new("yuima.qmleLevy.incr",Incr.Lev=Incr.Lev,
+ minusloglLevy = minusloglik,logL.Incr=-esti$value,
+ Data = yuima at data, yuima=res, Levydetails=esti)
+ vcovLevy <- function(myres, HessianEta_divTn, Gammahat0, sq=TRUE){
+ #myres <- res.VG2
+ DeltaX <- apply(myres at Data@original.data,2,diff)
+ myData<- myres at Data@original.data
+ CmatExpr<- myres at model@jump.coeff
+ ncolC <-length(myres at model@jump.coeff[[1]])
+ param <- myres at coef[myres at model@parameter at jump]
+ aexpr<- myres at model@drift
+ namedrift<-myres at model@parameter at drift
+ pardrif <- myres at coef[namedrift]
+ avect_exp<-myres at model@drift
+
+ Jac_Drift <- function(aexpr,namedrift,nobs=length(aexpr)){
+ lapply(X=c(1:nobs),FUN=function(X,namedrift,aexpr){
+ return(deriv(aexpr[X],namedrift))
+ },namedrift=namedrift,aexpr=aexpr)
+ }
+
+ Jac_DriftExp <- Jac_Drift(aexpr=avect_exp,namedrift)
+
+ FUNDum<-function(foo,myenv) sapply(foo, function(x,env) eval(x,envir = env), env= myenv)
+
+ h<-diff(time(myres at Data@zoo.data[[1]]))[1]
+
+ dummyf1n<- function(DeltaX, Cmat, h){
+ C2<-t(Cmat)%*%Cmat
+ dec <- chol(C2) # Eventualy Add a tryCatch
+ tmp <- t(DeltaX)%*%solve(C2)%*%DeltaX
+ logretval <- -h*sum(log(diag(dec))) - 0.5 * as.numeric(tmp)
+ return(logretval)
+ }
+ dummyInc<- function(DeltaX,Cmat, avect,h,sq=TRUE){
+ if(sq){
+ Incr <- solve(Cmat)%*%(DeltaX-avect*h)
+ }else{
+ C2<-t(Cmat)%*%Cmat
+ Incr <- solve(C2)%*%Cmat%*%(DeltaX-avect*h)
+ }
+ return(Incr)
+ }
+
+ dummyf2n<- function(DeltaX, avect, Cmat, h){
+ C2 <- t(Cmat)%*%Cmat
+ Incr <- DeltaX-h*avect
+ tmp <- t(Incr)%*%solve(C2)%*%Incr
+ logretval <- - 0.5/h * as.numeric(tmp)
+ return(logretval)
+ }
+ dumGrad_f2n<-function(DeltaX, avect, Jac_a,Cmat, h){
+ C2 <- t(Cmat)%*%Cmat
+ Incr <- DeltaX-h*avect
+ Grad<-t(Jac_a)%*%solve(C2)%*%Incr
+ return(Grad)
+ }
+
+ f1n_j<-function(x, myObsDelta, CmatExpr,h, myData, myres){
+ names(myData)<-myres at model@solve.variable
+ newenvJumpCoef <- list2env(as.list(c(x,myData)))
+ Cmat<-sapply(X=CmatExpr, FUN=FUNDum, myenv=newenvJumpCoef )
+ return(dummyf1n(DeltaX=myObsDelta, Cmat=Cmat, h=h))
+ }
+
+ Incr_Func <- function(x, myObsDelta, CmatExpr,avect_exp, h, myData, myres, sq=TRUE){
+ names(myData)<-myres at model@solve.variable
+ newenv <- list2env(as.list(c(x,myData)))
+ Cmat<-sapply(X=CmatExpr, FUN=FUNDum, myenv=newenv )
+ dd<-length(avect_exp)
+ avect<-numeric(length=dd)
+ for(j in c(1:dd)){
+ avect[j]<-eval(avect_exp[j],envir=newenv)
+ }
+ return(dummyInc(DeltaX= myObsDelta,Cmat, avect,h,sq=sq))
+ }
+
+
+ f2n_j <- function(x, parDiff, myObsDelta, CmatExpr, avect_exp, h, myData, myres){
+ names(myData)<-myres at model@solve.variable
+ newenvJumpCoef <- list2env(as.list(c(parDiff,myData)))
+ Cmat<-sapply(X=CmatExpr, FUN=FUNDum, myenv=newenvJumpCoef)
+ newenvDriftCoef <- list2env(as.list(c(x,myData)))
+ dd<-length(avect_exp)
+ avect<-numeric(length=dd)
+ for(j in c(1:dd)){
+ avect[j]<-eval(avect_exp[j],envir=newenvDriftCoef)
+ }
+ return(dummyf2n(DeltaX=myObsDelta, avect=avect, Cmat=Cmat, h=h))
+ }
+
+ Gradf2n_j <-function(x, parDiff, myObsDelta, CmatExpr,
+ avect_exp, Jac_DriftExp, h, myData, myres){
+ names(myData)<-myres at model@solve.variable
+ newenvJumpCoef <- list2env(as.list(c(parDiff,myData)))
+ Cmat<-sapply(X=CmatExpr, FUN=FUNDum, myenv=newenvJumpCoef)
+ newenvDriftCoef <- list2env(as.list(c(x,myData)))
+ dd<-length(avect_exp)
+ avect<-numeric(length=dd)
+ Jac_a<- matrix(0,dd,length(x))
+ for(j in c(1:dd)){
+ avect[j]<-eval(avect_exp[j],envir=newenvDriftCoef)
+ Jac_a[j,]<-attr(eval(Jac_DriftExp[[j]],envir=newenvDriftCoef),"gradient")
+ }
+ return(dumGrad_f2n(DeltaX=myObsDelta, avect=avect, Jac_a=Jac_a, Cmat=Cmat, h=h))
+ }
+
+ # f1n_j(x=param, myObsDelta=DeltaX[1,], CmatExpr=CmatExpr, h=h)
+ #i=1
+ # debug(dummyInc)
+ # Incr_Func(x=c(param,pardrif), myObsDelta=DeltaX[i,], CmatExpr,avect_exp, h, myData=myData[i,], myres, sq=TRUE)
+ # f1n_j(x=param, myObsDelta=DeltaX[i,], CmatExpr=CmatExpr, h=h, myData=myData[i,], myres=myres)
+ # f2n_j(x, parDiff=param, myObsDelta=DeltaX[i,], CmatExpr=CmatExpr, avect_exp=avect_exp,
+ # h=h, myData=myData[i,], myres=myres)
+ # Gradf2n_j(x, parDiff=param, myObsDelta=DeltaX[i,], CmatExpr,
+ # avect_exp, Jac_DriftExp, h, myData[i,], myres)
+ del <- 10^-3
+ dummy <- t(rep(1,length(param)))
+ myeta<- as.matrix(param)%*%dummy
+ dummyG <- t(rep(1,length(c(param,pardrif))))
+ globalmyeta <- as.matrix(c(param,pardrif))%*%dummyG
+ Incr <- del*diag(rep(1,dim(myeta)[1]))
+ GlobIncr <- del*diag(rep(1,dim(globalmyeta)[1]))
+ Sigma_gamma <- matrix(0 ,length(param),length(param))
+ Sigma_alpha <- matrix(0 ,length(pardrif),length(pardrif))
+ Sigma_algam <- matrix(0 ,length(pardrif),length(param))
+ histDeltaf1n_j_delta <-matrix(0 , dim(DeltaX)[1],length(param))
+ histDeltaf2n_j_delta <- matrix(0 , dim(DeltaX)[1],length(pardrif))
+ histb_incr <- array(0, c(dim(DeltaX)[2],dim(DeltaX)[1],length(c(param,pardrif))))
+ for(i in c(1:dim(DeltaX)[1])){
+ Deltaf1n_j_delta<-sapply(X=1:dim(myeta)[1],
+ FUN = function(X,myObsDelta, CmatExpr, h, myData, myres,del){
+ par1<-myeta[,X]+Incr[,X]
+ #par[oldyuima at model@measure$df at time.var]<-1
+ f1<-f1n_j(x=par1, myObsDelta, CmatExpr, h, myData, myres)
+ par2<-myeta[,X]-Incr[,X]
+ f2<-f1n_j(x=par2, myObsDelta, CmatExpr, h, myData, myres)
+ return((f1-f2)/(2*del))
+ },
+ myObsDelta=DeltaX[i,], CmatExpr=CmatExpr, h=h, myData=myData[i,],
+ myres=myres, del=del)
+
+ DeltaInc <- sapply(X=1:dim(GlobIncr)[1],
+ FUN = function(X, myObsDelta, CmatExpr,avect_exp, h,
+ myData, myres, del, sq){
+ par1<-globalmyeta[,X]+GlobIncr[,X]
+
+ # f1<-f1n_j(x=par1, myObsDelta, CmatExpr, h, myData, myres)
+ f1<-Incr_Func(x=par1, myObsDelta, CmatExpr,avect_exp, h, myData, myres, sq)
+ par2<-globalmyeta[,X]-GlobIncr[,X]
+ f2<-Incr_Func(x=par2, myObsDelta, CmatExpr,avect_exp, h, myData, myres, sq)
+ #f2<-f1n_j(x=par2, myObsDelta, CmatExpr, h, myData, myres)
+ return((f1-f2)/(2*del))
+ },
+ myObsDelta=DeltaX[i,], CmatExpr,avect_exp, h, myData=myData[i,], myres, del=del,sq=sq
+ )
+
+ histDeltaf1n_j_delta[i, ]<-Deltaf1n_j_delta
+
+ histb_incr[ , i, ] <- DeltaInc
+ Sigma_gamma <- Sigma_gamma + as.matrix(Deltaf1n_j_delta)%*%t(Deltaf1n_j_delta)
+
+ Deltaf2n_j_delta<- Gradf2n_j(x=pardrif, parDiff=param, myObsDelta=DeltaX[i,], CmatExpr,
+ avect_exp, Jac_DriftExp, h, myData[i,], myres)
+ histDeltaf2n_j_delta[i,]<-Deltaf2n_j_delta
+
+ Sigma_alpha <- Sigma_alpha + Deltaf2n_j_delta%*%t(Deltaf2n_j_delta)
+ Sigma_algam <- Sigma_algam + Deltaf2n_j_delta%*%t(Deltaf1n_j_delta)
+ #cat("\n",i)
+ }
+ Tn<-tail(time(myres at Data@original.data),1L)
+ Sigma_gamma0<-Sigma_gamma/Tn
+ Sigma_alpha0 <- Sigma_alpha/Tn
+ Sigma_algam0 <- Sigma_algam/Tn
+ Sigma0 <- cbind(rbind(Sigma_gamma0,Sigma_algam0),rbind(t(Sigma_algam0),Sigma_alpha0))
+
+ myLogDens <- function(x, myObsDelta, CmatExpr,h,myData, myres, pos){
+ return(f1n_j(x, myObsDelta[pos, ], CmatExpr,h,myData=myData[pos,], myres=myres))
+ }
+
+ myLogDens2<- function(x, parDiff, myObsDelta, CmatExpr, avect_exp, h, myData, myres, pos){
+ return(f2n_j(x, parDiff, myObsDelta[pos, ], CmatExpr, avect_exp, h, myData[pos,], myres))
+ #return(f1n_j(x, myObsDelta[pos, ], CmatExpr,h,myData=myData[pos,], myres=myres))
+ }
+
+
+ VectLogDens<-Vectorize(FUN=myLogDens,vectorize.args = "pos")
+ #VectLogDens2<-Vectorize(FUN=myLogDens2,vectorize.args = "pos")
+
+ mylogLik<- function(par,DeltaX,CmatExpr, myData, myres, h, Tn, nobs=dim(DeltaX)[1]){
+ term<-VectLogDens(x=par,
+ myObsDelta=DeltaX, CmatExpr=CmatExpr,h=h, myData=myData, myres=myres,pos=c(1:nobs))
+ return(sum(term)/Tn)
+ }
+
+ mylogLik2<- function(par, parDiff, myObsDelta, CmatExpr, avect_exp, h, Tn, myData, myres, nobs=dim(DeltaX)[1]){
+ # term<-VectLogDens2(x=par, parDiff, myObsDelta, CmatExpr, avect_exp, h, myData, myres,
+ # pos=c(1:nobs))
+ term<-0
+ for(j in c(1:nobs)){
+ term<-term+f2n_j(x=par, parDiff, myObsDelta[j, ], CmatExpr, avect_exp, h, myData[j,], myres)
+ }
+ return(term/Tn)
+ }
+
+ GradlogLik2<-function(par, parDiff, myObsDelta, CmatExpr, avect_exp, h, Tn, myData, myres, nobs=dim(DeltaX)[1]){
+ term<-matrix(0, length(par),1)
+ for(j in c(1:nobs)){
+ term <- term + Gradf2n_j(x=par, parDiff, myObsDelta=myObsDelta[j, ], CmatExpr,
+ avect_exp, Jac_DriftExp, h, myData=myData[j,], myres)
+ }
+ return(as.numeric(term)/Tn)
+ }
+
+ # Gammahat <- optimHess(par=param,fn=mylogLik,
+ # DeltaX=DeltaX,CmatExpr=CmatExpr,
+ # h=h, Tn=Tn, nobs=dim(DeltaX)[1],
+ # myData=myData, myres=myres)
+
+ # GammaAlfahat <- optimHess(par=pardrif,fn=mylogLik2,gr=GradlogLik2,
+ # parDiff=param, myObsDelta=DeltaX,
+ # CmatExpr=CmatExpr, avect_exp=avect_exp, h=h,
+ # Tn=Tn, nobs=dim(DeltaX)[1], myData=myData, myres=myres)
+ #
+ # Gammahat0 <- cbind(rbind(Gammahat,matrix(0,dim(GammaAlfahat)[1],dim(Gammahat)[2])),
+ # rbind(matrix(0,dim(Gammahat)[2],dim(GammaAlfahat)[1]),GammaAlfahat))
+
+ uhistDeltaf1n_j<- matrix(0,floor(Tn),length(param))
+ uhistDeltaf2n_j<- matrix(0,floor(Tn),length(pardrif))
+ aaa<-length(c(param,pardrif))
+ uhistb_incr <- array(0, c(dim(DeltaX)[2],floor(Tn),aaa))
+ for(l in 1:floor(Tn)){
+ #ures[l] <- sum(resi[(floor((l-1)/h)):(floor(l/h)-1)])
+ uhistDeltaf1n_j[l, ] <- colSums(histDeltaf1n_j_delta[(floor((l-1)/h)):(floor(l/h)-1),])
+ uhistDeltaf2n_j[l, ] <- colSums(histDeltaf2n_j_delta[(floor((l-1)/h)):(floor(l/h)-1),])
+ for(i in c(1:aaa)){
+ dumIn <- histb_incr[,(floor((l-1)/h)):(floor(l/h)-1),i]
+ uhistb_incr[,l,i]<- rowSums(dumIn)
+ }
+ }
+
+ nMeaspar<-myres at model@parameter at measure
+ mypar_meas0 <- myres at coef[nMeaspar]#res at coef[oldyuima at model@parameter at measure]
+ mypar_meas <- c(mypar_meas0,1)
+ ures<-myres at Incr.Lev@original.data
+ names(mypar_meas)<-c(nMeaspar,myres at model@measure$df at time.var)
+ fdataeta<-dens(object=myres at model@measure$df,x=ures,param=mypar_meas)# f(eps, eta)
+ del <- 10^-3
+ dummy <- t(rep(1,length(mypar_meas[myres at model@measure$df at param.measure])))
+ myeta<- as.matrix(mypar_meas[myres at model@measure$df at param.measure])%*%dummy
+ myetapert <- myeta+del*diag(rep(1,dim(myeta)[1]))
+ fdataetadelta<-sapply(X=1:dim(myeta)[1],FUN = function(X){
+ par<-myetapert[,X]
+ par[myres at model@measure$df at time.var]<-1
+ dens(object=myres at model@measure$df,x=ures,param=par)
+ }
+ )# f(eps, eta+delta)
+
+ term1<-1/(fdataeta)
+
+ dummyvecdel<-numeric(length=dim(ures)[2])
+ fdatadeltaeta<- matrix(0,floor(Tn),dim(ures)[2])
+ mixpartial <-array(0,c(length(nMeaspar),floor(Tn),dim(ures)[2]))
+ for(j in c(1:dim(ures)[2])){
+ dummyvecdel[j]<-del
+ fdatadeltaeta[,j]<- dens(object=myres at model@measure$df,x=ures+rep(1,dim(ures)[1])%*%t(dummyvecdel),param=mypar_meas)
+
+ fdatadeltaetadelta<-sapply(X=1:dim(myeta)[1],FUN = function(X){
+ par<-myetapert[,X]
+ par[myres at model@measure$df at time.var]<-1
+ dens(object=myres at model@measure$df,x=ures+rep(1,dim(ures)[1])%*%t(dummyvecdel),param=par)
+ }
+ ) # f(eps +deta, eta+delta)
+ dummyvecdel<-numeric(length=dim(ures)[2])
+ term2 <- fdatadeltaeta[,j]*term1%*%dummy
+ term2 <- fdatadeltaetadelta-term2*fdataetadelta
+ mixpartial[,,j]<-t((as.matrix(term1)%*%dummy)/del^2*term2/Tn)
+ }
+
+ DerMeta <- 1/del*(fdataetadelta - as.matrix(fdataeta)%*%rep(1,dim(fdataetadelta)[2]))*(as.matrix(term1)%*%rep(1,dim(fdataetadelta)[2]))
+ SigmaEta <- t(DerMeta)%*%DerMeta/Tn
+
+ minusloglik <- function(para){
+ para[length(para)+1]<-1
+ names(para)[length(para)]<-myres at model@time.variable
+ -sum(dens(object=myres at model@measure$df, x=ures, param = para, log = TRUE),
+ na.rm = T)/Tn
+ }
+ # HessianEta_divTn <- optimHess(par=mypar_meas0, fn=minusloglik)
+
+ Gammaeta_theta <- matrix(0,length(mypar_meas0), aaa)
+ for(t in c(1:floor(Tn))){
+ Gammaeta_theta<-Gammaeta_theta+mixpartial[,t,]%*%uhistb_incr[,t,]
+ }
+ GammaEta_Theta<- -1/Tn*Gammaeta_theta
+
+ Sigma_eta_theta<-t(DerMeta)%*%cbind(uhistDeltaf1n_j,uhistDeltaf2n_j)/Tn
+
+ Sigma <- cbind(rbind(Sigma0,Sigma_eta_theta),rbind(t(Sigma_eta_theta),SigmaEta))
+
+ GammaHAT<-cbind(rbind(Gammahat0, GammaEta_Theta), rbind(t(GammaEta_Theta), HessianEta_divTn))
+ InvGammaHAT <- solve(GammaHAT)
+ vcov <- InvGammaHAT %*% Sigma %*% InvGammaHAT/Tn
+
+ myres at vcov<- vcov
+ return(myres)
+ }
+
+ sq<-TRUE
+ result <- vcovLevy(result,HessianEta_divTn,myGamhat,sq=sq)
+ # colnames(result at vcov)<-names(res at fullcoef)
+ # rownames(result at vcov)<-names(res at fullcoef)
}
- colnames(res at vcov)<-names(res at fullcoef)
- #res at vcov<-rbind(res at vcov,matrix(NA,nrow=length(esti$par),ncol=dim(res at vcov)[2]))
- rownames(res at vcov)<-names(res at fullcoef)
- res at min<-c(res at min,esti$value)
- res at nobs<-c(res at nobs,length(Incr.Lev at zoo.data[[1]]))
+ # colnames(res at vcov)<-names(res at fullcoef)
+ # #res at vcov<-rbind(res at vcov,matrix(NA,nrow=length(esti$par),ncol=dim(res at vcov)[2]))
+ # rownames(res at vcov)<-names(res at fullcoef)
+ # res at min<-c(res at min,esti$value)
+ # res at nobs<-c(res at nobs,length(Incr.Lev at zoo.data[[1]]))
- result<- new("yuima.qmleLevy.incr",Incr.Lev=Incr.Lev,
- minusloglLevy = minusloglik,logL.Incr=-esti$value,
- Data = yuima at data, yuima=res, Levydetails=esti)
+
return(result)
}else{
dist <- substr(as.character(orig.mylaw$df$expr), 2, 10^3)
More information about the Yuima-commits
mailing list