[Yuima-commits] r278 - in pkg/yuima: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Mar 16 22:02:07 CET 2014
Author: lorenzo
Date: 2014-03-16 22:02:07 +0100 (Sun, 16 Mar 2014)
New Revision: 278
Modified:
pkg/yuima/R/AllClasses.R
pkg/yuima/R/qmle.R
pkg/yuima/man/qmle.Rd
Log:
Modified qmle for levy carma
Modified: pkg/yuima/R/AllClasses.R
===================================================================
--- pkg/yuima/R/AllClasses.R 2014-02-16 14:22:04 UTC (rev 277)
+++ pkg/yuima/R/AllClasses.R 2014-03-16 21:02:07 UTC (rev 278)
@@ -122,7 +122,7 @@
)
# Class yuima.carma.qmle
-setClass("yuima.carma.qmle",representation(Incr.Lev = "matrix",
+setClass("yuima.carma.qmle",representation(Incr.Lev = "ANY",
model = "yuima.carma"
),
contains="mle"
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2014-02-16 14:22:04 UTC (rev 277)
+++ pkg/yuima/R/qmle.R 2014-03-16 21:02:07 UTC (rev 278)
@@ -73,7 +73,7 @@
qmle <- function(yuima, start, method="BFGS", fixed = list(), print=FALSE,
- lower, upper, joint=FALSE, Est.Incr="Carma.IncPar", ...){
+ lower, upper, joint=FALSE, Est.Incr="Carma.IncPar",aggregation=TRUE, ...){
if(is(yuima at model, "yuima.carma")){
NoNeg.Noise<-FALSE
cat("\nStarting qmle for carma ... \n")
@@ -343,6 +343,7 @@
# env$X<-as.matrix(env$X[,1])
# assign("deltaX", matrix(0, n-1, d.size)[,1], envir=env)
env$X<-as.matrix(env$X[,1])
+# env$X<-na.omit(as.matrix(env$X[,1]))
env$deltaX<-as.matrix(env$deltaX[,1])
assign("time.obs",length(env$X),envir=env)
# env$time.obs<-length(env$X)
@@ -745,9 +746,10 @@
}
# INSERT HERE THE NECESSARY STEPS FOR FINDING THE PARAMETERS OF LEVY
if(Est.Incr=="Carma.Inc"){
+ inc.levy.fin<-zoo(inc.levy,tt,frequency=1/env$h)
carma_final_res<-new("yuima.carma.qmle", call = call, coef = coef, fullcoef = unlist(mycoef),
vcov = vcov, min = min, details = oout, minuslogl = minusquasilogl,
- method = method, Incr.Lev = inc.levy,
+ method = method, Incr.Lev = inc.levy.fin,
model = yuima at model)
return(carma_final_res)
}
@@ -804,29 +806,67 @@
return(carma_final_res)
}
- inc.levy1<-diff(cumsum(inc.levy)[seq(from=1,
- to=yuima at sampling@n[1],
- by=(yuima at sampling@n/yuima at sampling@Terminal)[1]
- )])
+ if(aggregation==TRUE){
+ inc.levy1<-diff(cumsum(inc.levy)[seq(from=1,
+ to=yuima at sampling@n[1],
+ by=(yuima at sampling@n/yuima at sampling@Terminal)[1]
+ )])
+ }else{
+ inc.levy1<-inc.levy
+ }
+
names.measpar<-c(name.int.dummy, names.measpar)
if(measurefunc=="dnorm"){
- result.Lev<-yuima.Estimation.CPN(Increment.lev=inc.levy1,param0=coef[ names.measpar],
+
+# result.Lev<-yuima.Estimation.CPN(Increment.lev=inc.levy1,param0=coef[ names.measpar],
+# fixed.carma=fixed.carma,
+# lower.carma=lower.carma,
+# upper.carma=upper.carma)
+
+ result.Lev<-yuima.Estimation.Lev(Increment.lev=inc.levy1,
+ param0=coef[ names.measpar],
fixed.carma=fixed.carma,
lower.carma=lower.carma,
- upper.carma=upper.carma)
+ upper.carma=upper.carma,
+ measure=measurefunc,
+ measure.type=model at measure.type,
+ dt=env$h,
+ aggregation=aggregation)
+
}
if(measurefunc=="dgamma"){
- result.Lev<-yuima.Estimation.CPGam(Increment.lev=inc.levy1,param0=coef[ names.measpar],
- fixed.carma=fixed.carma,
- lower.carma=lower.carma,
- upper.carma=upper.carma)
+# result.Lev<-yuima.Estimation.CPGam(Increment.lev=inc.levy1,param0=coef[ names.measpar],
+# fixed.carma=fixed.carma,
+# lower.carma=lower.carma,
+# upper.carma=upper.carma)
+
+ result.Lev<-yuima.Estimation.Lev(Increment.lev=inc.levy1,
+ param0=coef[ names.measpar],
+ fixed.carma=fixed.carma,
+ lower.carma=lower.carma,
+ upper.carma=upper.carma,
+ measure=measurefunc,
+ measure.type=model at measure.type,
+ dt=env$h,
+ aggregation=aggregation)
}
if(measurefunc=="dexp"){
- result.Lev<-yuima.Estimation.CPExp(Increment.lev=inc.levy1,param0=coef[ names.measpar],
- fixed.carma=fixed.carma,
- lower.carma=lower.carma,
- upper.carma=upper.carma)
+# result.Lev<-yuima.Estimation.CPExp(Increment.lev=inc.levy1,param0=coef[ names.measpar],
+# fixed.carma=fixed.carma,
+# lower.carma=lower.carma,
+# upper.carma=upper.carma)
+
+ result.Lev<-yuima.Estimation.Lev(Increment.lev=inc.levy1,
+ param0=coef[ names.measpar],
+ fixed.carma=fixed.carma,
+ lower.carma=lower.carma,
+ upper.carma=upper.carma,
+ measure=measurefunc,
+ measure.type=model at measure.type,
+ dt=env$h,
+ aggregation=aggregation)
+
}
Inc.Parm<-result.Lev$estLevpar
IncVCOV<-result.Lev$covLev
@@ -867,13 +907,15 @@
model = yuima at model)
return(carma_final_res)
}
-
- inc.levy1<-diff(cumsum(inc.levy)[seq(from=1,
- to=yuima at sampling@n[1],
- by=(yuima at sampling@n/yuima at sampling@Terminal)[1]
- )])
-
-
+ if(aggregation==TRUE){
+ inc.levy1<-diff(cumsum(inc.levy)[seq(from=1,
+ to=yuima at sampling@n[1],
+ by=(yuima at sampling@n/yuima at sampling@Terminal)[1]
+ )])
+ }else{
+ inc.levy1<-inc.levy
+ }
+
if(measurefunc=="rIG"){
# result.Lev<-list(estLevpar=coef[ names.measpar],
@@ -881,20 +923,39 @@
# length(coef[ names.measpar]),
# length(coef[ names.measpar]))
# )
- result.Lev<-yuima.Estimation.IG(Increment.lev=inc.levy1,param0=coef[ names.measpar],
- fixed.carma=fixed.carma,
- lower.carma=lower.carma,
- upper.carma=upper.carma)
+# result.Lev<-yuima.Estimation.IG(Increment.lev=inc.levy1,param0=coef[ names.measpar],
+# fixed.carma=fixed.carma,
+# lower.carma=lower.carma,
+# upper.carma=upper.carma)
- # result.Levy<-gigFit(inc.levy)
+ result.Lev<-yuima.Estimation.Lev(Increment.lev=inc.levy1,
+ param0=coef[ names.measpar],
+ fixed.carma=fixed.carma,
+ lower.carma=lower.carma,
+ upper.carma=upper.carma,
+ measure=measurefunc,
+ measure.type=model at measure.type,
+ dt=env$h,
+ aggregation=aggregation)
+ # result.Levy<-gigFit(inc.levy)
# Inc.Parm<-coef(result.Levy)
# IncVCOV<--solve(gigHessian(inc.levy, param=Inc.Parm))
}
if(measurefunc=="rNIG"){
- result.Lev<-yuima.Estimation.NIG(Increment.lev=inc.levy1,param0=coef[ names.measpar],
- fixed.carma=fixed.carma,
- lower.carma=lower.carma,
- upper.carma=upper.carma)
+# result.Lev<-yuima.Estimation.NIG(Increment.lev=inc.levy1,param0=coef[ names.measpar],
+# fixed.carma=fixed.carma,
+# lower.carma=lower.carma,
+# upper.carma=upper.carma)
+
+ result.Lev<-yuima.Estimation.Lev(Increment.lev=inc.levy1,
+ param0=coef[ names.measpar],
+ fixed.carma=fixed.carma,
+ lower.carma=lower.carma,
+ upper.carma=upper.carma,
+ measure=measurefunc,
+ measure.type=model at measure.type,
+ dt=env$h,
+ aggregation=aggregation)
}
if(measurefunc=="rbgamma"){
result.Lev<-list(estLevpar=coef[ names.measpar],
@@ -904,10 +965,21 @@
)
}
if(measurefunc=="rngamma"){
- result.Lev<-yuima.Estimation.VG(Increment.lev=inc.levy1,param0=coef[ names.measpar],
- fixed.carma=fixed.carma,
- lower.carma=lower.carma,
- upper.carma=upper.carma)
+# result.Lev<-yuima.Estimation.VG(Increment.lev=inc.levy1,param0=coef[ names.measpar],
+# fixed.carma=fixed.carma,
+# lower.carma=lower.carma,
+# upper.carma=upper.carma)
+
+ result.Lev<-yuima.Estimation.Lev(Increment.lev=inc.levy1,
+ param0=coef[ names.measpar],
+ fixed.carma=fixed.carma,
+ lower.carma=lower.carma,
+ upper.carma=upper.carma,
+ measure=measurefunc,
+ measure.type=model at measure.type,
+ dt=env$h,
+ aggregation=aggregation)
+
}
Inc.Parm<-result.Lev$estLevpar
@@ -951,13 +1023,14 @@
# carma_final_res<-list(mle=final_res,Incr=inc.levy,model=yuima)
if(Est.Incr=="Carma.IncPar"){
- carma_final_res<-new("yuima.carma.qmle", call = call, coef = coef, fullcoef = unlist(mycoef),
+ inc.levy.fin<-zoo(inc.levy,tt,frequency=1/env$h)
+ carma_final_res<-new("yuima.carma.qmle", call = call, coef = coef, fullcoef = unlist(coef),
vcov = cov, min = min, details = oout, minuslogl = minusquasilogl,
- method = method, Incr.Lev = inc.levy,
+ method = method, Incr.Lev = inc.levy.fin,
model = yuima at model)
}else{
- if(Est.Incr=="Carma.IncPar"){
- carma_final_res<-new("mle", call = call, coef = coef, fullcoef = unlist(mycoef),
+ if(Est.Incr=="Carma.Par"){
+ carma_final_res<-new("mle", call = call, coef = coef, fullcoef = unlist(coef),
vcov = cov, min = min, details = oout, minuslogl = minusquasilogl,
method = method)
}
@@ -989,8 +1062,8 @@
# assign("deltaX", matrix(0, n-1, d.size)[,1], envir=env)
env$X<-as.matrix(env$X[,1])
env$deltaX<-as.matrix(env$deltaX[,1])
- assign("time.obs",length(env$X),envir=env)
- #env$time.obs<-length(env$X)
+ #assign("time.obs",length(env$X),envir=env)
+ env$time.obs<-length(env$X)
#p <-yuima at model@info at p
assign("p", yuima at model@info at p, envir=env)
assign("q", yuima at model@info at q, envir=env)
@@ -1397,16 +1470,25 @@
ATrans<-t(A)
matrixV<-matrix(0,p,p)
+ #l.dummy<-c(rep(0,p-1),1)
l<-rbind(matrix(rep(0,p-1),p-1,1),1)
+ #l<-matrix(l.dummy,p,1)
+ #lTrans<-matrix(l.dummy,1,p)
lTrans<-t(l)
-
- elForVInf<-list(A=A,
- ATrans=ATrans,
- lTrans=lTrans,
- l=l,
- matrixV=matrixV,
- sigma=sigma)
-
+ elForVInf<-new.env()
+ elForVInf$A<-A
+ elForVInf$ATrans<-ATrans
+ elForVInf$lTrans<-lTrans
+ elForVInf$l<-l
+ elForVInf$matrixV<-matrixV
+ elForVInf$sigma<-sigma
+# elForVInf<-list(A=A,
+# ATrans=ATrans,
+# lTrans=lTrans,
+# l=l,
+# matrixV=matrixV,
+# sigma=sigma)
+#
V_inf_vect<-nlm(yuima.Vinfinity, v,
elForVInf = elForVInf)$estimate
# V_inf_vect<-nlminb(start=v,objective=yuima.Vinfinity, elForVInf = elForVInf)$par
@@ -1427,7 +1509,7 @@
# set
#statevar<-statevar0
- # SigMatr<-expA%*%V_inf%*%t(expA)+Qmatr
+ # SigMatr<-expA%*%V_inf%*%t(expA)+Qmatr
#SigMatr<-Qmatr
#SigMatr<-V_inf
@@ -1436,6 +1518,7 @@
loglstar <- 0
loglstar1 <- 0
+# zcT<-matrix(bvector,p,1)
zcT<-t(zc)
for(t in 1:times.obs){
# prediction
@@ -1443,20 +1526,16 @@
SigMatr<-expA%*%SigMatr%*%expAT+Qmatr
# forecast
Uobs<-y[t]-zc%*%statevar
- sd_2<-zc%*%SigMatr%*%zcT
+ dum.zc<-zc%*%SigMatr
+ sd_2<-dum.zc%*%zcT
+ # sd_2<-zc%*%SigMatr%*%zcT
Inv_sd_2<-1/sd_2
#correction
- Kgain<-SigMatr%*%zcT%*%Inv_sd_2
- # SigMatr<-SigMatr-Kgain%*%as.numeric(sd_2)%*%t(Kgain)
-
+ Kgain<-SigMatr%*%zcT%*%Inv_sd_2
statevar<-statevar+Kgain%*%Uobs
- SigMatr<-SigMatr-Kgain%*%zc%*%SigMatr
- # construction of log-likelihood
- # loglstar1<-loglstar1+log(dnorm(as.numeric(Uobs), mean = 0, sd = sqrt(as.numeric(sd_2))))
- # sdsig<-sqrt(as.numeric(sd_2))
- # term_int<--0.5*(log((2*pi)*det(sd_2))+t(Uobs)%*%Inv_sd_2%*%Uobs)
- #term_int<--0.5*(log((2*pi)*sd_2)+t(Uobs)%*%Inv_sd_2%*%Uobs)
- term_int<--0.5*(log(sd_2)+Uobs^2*Inv_sd_2)
+ #SigMatr<-SigMatr-Kgain%*%zc%*%SigMatr
+ SigMatr<-SigMatr-Kgain%*%dum.zc
+ term_int<--0.5*(log(sd_2)+Uobs%*%Uobs%*%Inv_sd_2)
loglstar<-loglstar+term_int
}
return(list(loglstar=(loglstar-0.5*log(2*pi)*times.obs),s2hat=sd_2))
@@ -1912,696 +1991,378 @@
# Plot Method for yuima.carma.qmle
setMethod("plot",signature(x="yuima.carma.qmle"),
function(x, ...){
- plot(x at Incr.Lev, ...)
+ Time<-index(x at Incr.Lev)
+ Incr.L<-coredata(x at Incr.Lev)
+ if(is.complex(Incr.L)){
+ yuima.warn("Complex increments. We plot only the real part")
+ Incr.L<-Re(Incr.L)
+ }
+ plot(x=Time,y=Incr.L, ...)
}
)
# Utilities for estimation of levy in continuous arma model
-# Normal Inverse Gaussian
+#Density code for compound poisson
-yuima.Estimation.NIG<-function(Increment.lev,param0,
- fixed.carma=fixed.carma,
- lower.carma=lower.carma,
- upper.carma=upper.carma){
-
- minusloglik.dNIG<-function(par,data){
- alpha<-par[1]
- beta<-par[2]
- delta<-par[3]
- mu<-par[4]
- -sum(log(dNIG(data,alpha,beta,delta,mu)))
- }
-
- data<-Increment.lev
-
- # Only one problem
-
-
- ui<-rbind(c(1, -1, 0, 0),c(1, 1, 0, 0),c(1, 0, 0, 0),c(0, 0, 1, 0))
- ci<-c(0,0,0,10^(-6))
+#CPN
- if(!is.null(lower.carma)){
- lower.con<-matrix(0,length(lower.carma),length(param0))
- rownames(lower.con)<-names(lower.carma)
- colnames(lower.con)<-names(param0)
- numb.lower<-length(lower.carma)
- lower.con[names(lower.carma),names(lower.carma)]<-1*diag(numb.lower)
- dummy.lower.names<-paste0(names(lower.carma),".lower")
- rownames(lower.con)<-dummy.lower.names
- names(lower.carma)<-dummy.lower.names
- ui<-rbind(ui,lower.con)
- ci<-c(ci,lower.carma)
- #idx.lower.carma<-match(names(lower.carma),names(param0))
- }
- if(!is.null(upper.carma)){
- upper.con<-matrix(0,length(upper.carma),length(param0))
- rownames(upper.con)<-names(upper.carma)
- colnames(upper.con)<-names(param0)
- numb.upper<-length(upper.carma)
- upper.con[names(upper.carma),names(upper.carma)]<--1*diag(numb.upper)
- dummy.upper.names<-paste0(names(upper.carma),".upper")
- rownames(upper.con)<-dummy.upper.names
- names(upper.carma)<-dummy.upper.names
- ui<-rbind(ui,upper.con)
- ci<-c(ci,-upper.carma)
- }
- if(!is.null(fixed.carma)){
- names.fixed<-names(fixed.carma)
- numb.fixed<-length(fixed.carma)
- fixed.con<-matrix(0,length(fixed.carma),length(param0))
- rownames(fixed.con)<-names(fixed.carma)
- colnames(fixed.con)<-names(param0)
- fixed.con.bis<-fixed.con
- fixed.con[names(fixed.carma),names(fixed.carma)]<--1*diag(numb.fixed)
- fixed.con.bis[names(fixed.carma),names(fixed.carma)]<-1*diag(numb.fixed)
- dummy.fixed.names<-paste0(names(fixed.carma),".fixed.u")
- dummy.fixed.bis.names<-paste0(names(fixed.carma),".fixed.l")
- rownames(fixed.con)<-dummy.fixed.names
- rownames(fixed.con.bis)<-dummy.fixed.bis.names
- names(fixed.carma)<-dummy.fixed.names
- ui<-rbind(ui,fixed.con,fixed.con.bis)
- ci<-c(ci,-fixed.carma-10^-6,fixed.carma-10^-6)
- #ci<-c(ci,-fixed.carma,fixed.carma)
- }
-
-
- firs.prob<-tryCatch(constrOptim(theta=param0,
- f=minusloglik.dNIG,grad=NULL,ui=ui,ci=ci,data=data),
- error=function(theta){NULL})
-
- lengpar<-length(param0)
- paramLev<-NA*c(1:length(lengpar))
-
- if(!is.null(firs.prob)){
- paramLev<-firs.prob$par
- names(paramLev)<-names(param0)
- if(!is.null(fixed.carma)){
- paramLev[names.fixed]<-fixed.carma
- names(paramLev)<-names(param0)
+dCPN<-function(x,lambda,mu,sigma){
+ a<-min(mu-100*sigma,min(x)-1)
+ b<-max(mu+100*sigma,max(x)+1)
+ ChFunToDens.CPN <- function(n, a, b, lambda, mu, sigma) {
+ i <- 0:(n-1) # Indices
+ dx <- (b-a)/n # Step size, for the density
+ x <- a + i * dx # Grid, for the density
+ dt <- 2*pi / ( n * dx ) # Step size, frequency space
+ c <- -n/2 * dt # Evaluate the characteristic function on [c,d]
+ d <- n/2 * dt # (center the interval on zero)
+ t <- c + i * dt # Grid, frequency space
+ charact.CPN<-function(t,lambda,mu,sigma){
+ normal.y<-exp(1i*t*mu-sigma^2*t^2/2)
+ y<-exp(lambda*(normal.y-1))
}
- }else{warning("the start value for levy measure is outside of the admissible region")}
-
- NIG.hessian<-function (data,params){
- logLik.NIG <- function(params) {
-
- alpha<-params[1]
- beta<-params[2]
- delta<-params[3]
- mu<-params[4]
-
- return(sum(log(dNIG(data,alpha,beta,delta,mu))))
- }
- hessian<-optimHess(par=params, fn=logLik.NIG)
- cov<--solve(hessian)
- return(cov)
+ phi_t <- charact.CPN(t,lambda,mu,sigma)
+ X <- exp( -(0+1i) * i * dt * a ) * phi_t
+ Y <- fft(X)
+ density <- dt / (2*pi) * exp( - (0+1i) * c * x ) * Y
+ data.frame(
+ i = i,
+ t = t,
+ characteristic_function = phi_t,
+ x = x,
+ density = Re(density)
+ )
}
-
- if(is.na(paramLev)){
- covLev<-matrix(NA,length(paramLev),length(paramLev))
- }else{
- covLev<-NIG.hessian(data=as.numeric(data),params=paramLev)
- rownames(covLev)<-names(paramLev)
- if(!is.null(fixed.carma)){
- covLev[names.fixed,]<-matrix(NA,numb.fixed,lengpar)
+ invFFT<-ChFunToDens.CPN(lambda=lambda,mu=mu,sigma=sigma,n=2^10,a=a,b=b)
+ dens<-approx(invFFT$x,invFFT$density,x)
+ return(dens$y)
+}
+
+# CExp
+
+dCPExp<-function(x,lambda,rate){
+ a<-10^-6
+ b<-max(1/rate*10 +1/rate^2*10 ,max(x[!is.na(x)])+1)
+ ChFunToDens.CPExp <- function(n, a, b, lambda, rate) {
+ i <- 0:(n-1) # Indices
+ dx <- (b-a)/n # Step size, for the density
+ x <- a + i * dx # Grid, for the density
+ dt <- 2*pi / ( n * dx ) # Step size, frequency space
+ c <- -n/2 * dt # Evaluate the characteristic function on [c,d]
+ d <- n/2 * dt # (center the interval on zero)
+ t <- c + i * dt # Grid, frequency space
+ charact.CPExp<-function(t,lambda,rate){
+ normal.y<-(rate/(1-1i*t))
+ # exp(1i*t*mu-sigma^2*t^2/2)
+ y<-exp(lambda*(normal.y-1))
}
- colnames(covLev)<-names(paramLev)
- if(!is.null(fixed.carma)){
- covLev[,names.fixed]<-matrix(NA,lengpar,numb.fixed)
- }
+ phi_t <- charact.CPExp(t,lambda,rate)
+ X <- exp( -(0+1i) * i * dt * a ) * phi_t
+ Y <- fft(X)
+ density <- dt / (2*pi) * exp( - (0+1i) * c * x ) * Y
+ data.frame(
+ i = i,
+ t = t,
+ characteristic_function = phi_t,
+ x = x,
+ density = Re(density)
+ )
}
- results<-list(estLevpar=paramLev,covLev=covLev)
- return(results)
+ invFFT<-ChFunToDens.CPExp(lambda=lambda,rate=rate,n=2^10,a=a,b=b)
+ dens<-approx(invFFT$x[!is.na(invFFT$density)],invFFT$density[!is.na(invFFT$density)],x)
+ return(dens$y[!is.na(dens$y)])
}
+# CGamma
+dCPGam<-function(x,lambda,shape,scale){
+ a<-10^-6
+ b<-max(shape*scale*10 +shape*scale^2*10 ,max(x[!is.na(x)])+1)
+ ChFunToDens.CPGam <- function(n, a, b, lambda, shape,scale) {
+ i <- 0:(n-1) # Indices
+ dx <- (b-a)/n # Step size, for the density
+ x <- a + i * dx # Grid, for the density
+ dt <- 2*pi / ( n * dx ) # Step size, frequency space
+ c <- -n/2 * dt # Evaluate the characteristic function on [c,d]
+ d <- n/2 * dt # (center the interval on zero)
+ t <- c + i * dt # Grid, frequency space
+ charact.CPGam<-function(t,lambda,shape,scale){
+ normal.y<-(1-1i*t*scale)^(-shape)
+ # exp(1i*t*mu-sigma^2*t^2/2)
+ y<-exp(lambda*(normal.y-1))
+ }
+ phi_t <- charact.CPGam(t,lambda,shape,scale)
+ X <- exp( -(0+1i) * i * dt * a ) * phi_t
+ Y <- fft(X)
+ density <- dt / (2*pi) * exp( - (0+1i) * c * x ) * Y
+ data.frame(
+ i = i,
+ t = t,
+ characteristic_function = phi_t,
+ x = x,
+ density = Re(density)
+ )
+ }
+ invFFT<-ChFunToDens.CPGam(lambda=lambda,shape=shape,scale=scale,n=2^10,a=a,b=b)
+ dens<-approx(invFFT$x[!is.na(invFFT$density)],invFFT$density[!is.na(invFFT$density)],x)
+ return(dens$y[!is.na(dens$y)])
+}
-# Variance Gaussian
-yuima.Estimation.VG<-function(Increment.lev,param0,
- fixed.carma=fixed.carma,
- lower.carma=lower.carma,
- upper.carma=upper.carma){
-
- minusloglik.dVG<-function(par,data){
- lambda<-par[1]
- alpha<-par[2]
- beta<-par[3]
- mu<-par[4]
- -sum(log(dngamma(data,lambda,alpha,beta,mu)))
- }
-
- data<-Increment.lev
-
- ui<-rbind(c(1,0, 0, 0),c(0, 1, 1, 0),c(0, 1,-1, 0),c(0, 1,0, 0))
- ci<-c(10^-6,10^-6,10^(-6), 0)
-
- if(!is.null(lower.carma)){
- lower.con<-matrix(0,length(lower.carma),length(param0))
- rownames(lower.con)<-names(lower.carma)
- colnames(lower.con)<-names(param0)
- numb.lower<-length(lower.carma)
- lower.con[names(lower.carma),names(lower.carma)]<-1*diag(numb.lower)
- dummy.lower.names<-paste0(names(lower.carma),".lower")
- rownames(lower.con)<-dummy.lower.names
- names(lower.carma)<-dummy.lower.names
- ui<-rbind(ui,lower.con)
- ci<-c(ci,lower.carma)
- #idx.lower.carma<-match(names(lower.carma),names(param0))
- }
- if(!is.null(upper.carma)){
- upper.con<-matrix(0,length(upper.carma),length(param0))
- rownames(upper.con)<-names(upper.carma)
- colnames(upper.con)<-names(param0)
- numb.upper<-length(upper.carma)
- upper.con[names(upper.carma),names(upper.carma)]<--1*diag(numb.upper)
- dummy.upper.names<-paste0(names(upper.carma),".upper")
- rownames(upper.con)<-dummy.upper.names
- names(upper.carma)<-dummy.upper.names
- ui<-rbind(ui,upper.con)
- ci<-c(ci,-upper.carma)
- }
- if(!is.null(fixed.carma)){
- names.fixed<-names(fixed.carma)
- numb.fixed<-length(fixed.carma)
- fixed.con<-matrix(0,length(fixed.carma),length(param0))
- rownames(fixed.con)<-names(fixed.carma)
- colnames(fixed.con)<-names(param0)
- fixed.con.bis<-fixed.con
- fixed.con[names(fixed.carma),names(fixed.carma)]<--1*diag(numb.fixed)
- fixed.con.bis[names(fixed.carma),names(fixed.carma)]<-1*diag(numb.fixed)
- dummy.fixed.names<-paste0(names(fixed.carma),".fixed.u")
- dummy.fixed.bis.names<-paste0(names(fixed.carma),".fixed.l")
- rownames(fixed.con)<-dummy.fixed.names
- rownames(fixed.con.bis)<-dummy.fixed.bis.names
- names(fixed.carma)<-dummy.fixed.names
- ui<-rbind(ui,fixed.con,fixed.con.bis)
- ci<-c(ci,-fixed.carma-10^-6,fixed.carma-10^-6)
- #ci<-c(ci,-fixed.carma,fixed.carma)
- }
-
-
- firs.prob<-tryCatch(constrOptim(theta=param0,
- f=minusloglik.dVG,grad=NULL,ui=ui,ci=ci,data=data),
- error=function(theta){NULL})
-
- lengpar<-length(param0)
- paramLev<-NA*c(1:length(lengpar))
-
- if(!is.null(firs.prob)){
- paramLev<-firs.prob$par
- names(paramLev)<-names(param0)
- if(!is.null(fixed.carma)){
- paramLev[names.fixed]<-fixed.carma
- names(paramLev)<-names(param0)
+minusloglik.Lev<-function(par,env){
+ if(env$measure.type=="code"){
+ if(env$measure=="rNIG"){
+ alpha<-par[1]
+ beta<-par[2]
+ delta<-par[3]
+ mu<-par[4]
+ -sum(log(dNIG(env$data,alpha,beta,delta,mu)))
+ }else{
+ if(env$measure=="rngamma"){
+ lambda<-par[1]
+ alpha<-par[2]
+ beta<-par[3]
+ mu<-par[4]
+ -sum(log(dngamma(env$data,lambda,alpha,beta,mu)))
+ }else{
+ if(env$measure=="rIG"){
+ delta<-par[1]
+ gamma<-par[2]
+ f<-dIG(env$data,delta,gamma)
+ v<-log(as.numeric(na.omit(f)))
+ v1<-v[!is.infinite(v)]
+ -sum(v1)
+ }
+ }
}
- }
-
-
- VG.hessian<-function (data,params){
- logLik.VG <- function(params) {
-
- lambda<-params[1]
- alpha<-params[2]
- beta<-params[3]
- mu<-params[4]
-
- return(sum(log(dngamma(data,lambda,alpha,beta,mu))))
- }
- # hessian <- tsHessian(param = params, fun = logLik.VG)
- #hessian<-optimHess(par, fn, gr = NULL,data=data)
- hessian<-optimHess(par=params, fn=logLik.VG)
- cov<--solve(hessian)
- return(cov)
- }
-
- if(is.na(paramLev)){
- covLev<-matrix(NA,length(paramLev),length(paramLev))
}else{
- covLev<-VG.hessian(data=as.numeric(data),params=paramLev)
- rownames(covLev)<-names(paramLev)
- if(!is.null(fixed.carma)){
- covLev[names.fixed,]<-matrix(NA,numb.fixed,lengpar)
+ if(env$measure=="dnorm"){
+ lambda<-par[1]
+ mu<-par[2]
+ sigma<-par[3]
+ f<-dCPN(env$data,lambda,mu,sigma)
+ v<-log(as.numeric(na.omit(f)))
+ v1<-v[!is.infinite(v)]
+ -sum(v1)
+ }else{
+ if(env$measure=="dexp"){
+ lambda<-par[1]
+ rate<-par[2]
+ # -sum(log(dCPExp(env$data,lambda,rate)))
+
+ f<-dCPExp(env$data,lambda,rate)
+ v<-log(as.numeric(na.omit(f)))
+ v1<-v[!is.infinite(v)]
+ -sum(v1)
+
+ }else{
+ if(env$measure=="dgamma"){
+ lambda<-par[1]
+ shape<-par[2]
+ scale<-par[3]
+# -sum(log(dCPGam(env$data,lambda,shape,scale)))
+
+ f<-dCPGam(env$data,lambda,shape,scale)
+ v<-log(as.numeric(na.omit(f)))
+ v1<-v[!is.infinite(v)]
+ -sum(v1)
+
+
+ }
+ }
}
- colnames(covLev)<-names(paramLev)
- if(!is.null(fixed.carma)){
- covLev[,names.fixed]<-matrix(NA,lengpar,numb.fixed)
- }
}
- results<-list(estLevpar=paramLev,covLev=covLev)
- return(results)
}
-# Inverse Gaussian
-yuima.Estimation.IG<-function(Increment.lev,param0,
- fixed.carma=fixed.carma,
- lower.carma=lower.carma,
- upper.carma=upper.carma){
-
- minusloglik.dIG<-function(par,data){
- delta<-par[1]
- gamma<-par[2]
- -sum(log(dIG(data,delta,gamma)))
- }
-
- data<-Increment.lev
-
- ui<-rbind(c(1,0),c(0, 1))
- ci<-c(10^-6,10^-6)
-
- if(!is.null(lower.carma)){
- lower.con<-matrix(0,length(lower.carma),length(param0))
- rownames(lower.con)<-names(lower.carma)
- colnames(lower.con)<-names(param0)
- numb.lower<-length(lower.carma)
- lower.con[names(lower.carma),names(lower.carma)]<-1*diag(numb.lower)
- dummy.lower.names<-paste0(names(lower.carma),".lower")
- rownames(lower.con)<-dummy.lower.names
- names(lower.carma)<-dummy.lower.names
- ui<-rbind(ui,lower.con)
- ci<-c(ci,lower.carma)
- #idx.lower.carma<-match(names(lower.carma),names(param0))
- }
- if(!is.null(upper.carma)){
- upper.con<-matrix(0,length(upper.carma),length(param0))
- rownames(upper.con)<-names(upper.carma)
- colnames(upper.con)<-names(param0)
- numb.upper<-length(upper.carma)
- upper.con[names(upper.carma),names(upper.carma)]<--1*diag(numb.upper)
- dummy.upper.names<-paste0(names(upper.carma),".upper")
- rownames(upper.con)<-dummy.upper.names
- names(upper.carma)<-dummy.upper.names
- ui<-rbind(ui,upper.con)
- ci<-c(ci,-upper.carma)
- }
- if(!is.null(fixed.carma)){
- names.fixed<-names(fixed.carma)
- numb.fixed<-length(fixed.carma)
- fixed.con<-matrix(0,length(fixed.carma),length(param0))
- rownames(fixed.con)<-names(fixed.carma)
- colnames(fixed.con)<-names(param0)
- fixed.con.bis<-fixed.con
- fixed.con[names(fixed.carma),names(fixed.carma)]<--1*diag(numb.fixed)
- fixed.con.bis[names(fixed.carma),names(fixed.carma)]<-1*diag(numb.fixed)
- dummy.fixed.names<-paste0(names(fixed.carma),".fixed.u")
- dummy.fixed.bis.names<-paste0(names(fixed.carma),".fixed.l")
- rownames(fixed.con)<-dummy.fixed.names
- rownames(fixed.con.bis)<-dummy.fixed.bis.names
- names(fixed.carma)<-dummy.fixed.names
- ui<-rbind(ui,fixed.con,fixed.con.bis)
- ci<-c(ci,-fixed.carma-10^-6,fixed.carma-10^-6)
- #ci<-c(ci,-fixed.carma,fixed.carma)
- }
+Lev.hessian<-function (params,env){
+ logLik.Lev <- function(params){
+ if(env$measure.type=="code"){
+ if(env$measure=="rNIG"){
+ alpha<-params[1]
+ beta<-params[2]
+ delta<-params[3]
+ mu<-params[4]
+ # return(sum(log(dNIG(env$data,alpha,beta,delta,mu))))
+ f<-dNIG(env$data,alpha,beta,delta,mu)
+ v<-log(as.numeric(na.omit(f)))
+ v1<-v[!is.infinite(v)]
+ return(sum(v1))
+ }else{
+ if(env$measure=="rngamma"){
+ lambda<-params[1]
+ alpha<-params[2]
+ beta<-params[3]
+ mu<-params[4]
+ #return(sum(log(dngamma(env$data,lambda,alpha,beta,mu))))
+ f<-dngamma(env$data,lambda,alpha,beta,mu)
+ v<-log(as.numeric(na.omit(f)))
+ v1<-v[!is.infinite(v)]
+ return(sum(v1))
+ }else{
+ if(env$measure=="rIG"){
+ delta<-params[1]
+ gamma<-params[2]
+ f<-dIG(env$data,delta,gamma)
+ v<-log(as.numeric(na.omit(f)))
+ v1<-v[!is.infinite(v)]
+ return(sum(v1))
+ }else{
+ if(env$measure=="rgamma"){
- firs.prob<-tryCatch(constrOptim(theta=param0,
- f=minusloglik.dIG,
- grad=NULL,
- ui=ui,
- ci=ci,
- data=data),
- error=function(theta){NULL})
-
- lengpar<-length(param0)
- paramLev<-NA*c(1:length(lengpar))
- if(!is.null(firs.prob)){
- paramLev<-firs.prob$par
- names(paramLev)<-names(param0)
- if(!is.null(fixed.carma)){
- paramLev[names.fixed]<-fixed.carma
- names(paramLev)<-names(param0)
+ shape<-params[1]
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/yuima -r 278
More information about the Yuima-commits
mailing list