fracdiff.var(): did "FIXME"; testing "search" h; update tests

Author: mmaechler

.Rbuildignore

@@ -2,6 +2,8 @@ Archive
 Calling
 DESCRIPTION_
 Done
+M1mac
+00.*winb
 .*\.zip
 .*\.rar
 .*\.tgz
@@ -11,10 +13,11 @@ qed.*\.pdf
 codes_for_.*
 ^filters\.R
 simul-compare.*
-tests/.*linux
 tests/.*Rout-[36][24]b
-tests/windows.*
+tests/.*Rout.*_
 tests/ex-Vinod.*Rout
+tests/.*linux
+tests/windows.*
 src/.*\.c\+
 src/00-.*
 src/ftn-struc

ChangeLog

@@ -1,3 +1,9 @@
+2024-03-13  Martin Maechler  <[email protected]>
+
+	* DESCRIPTION (Version): 1.5-4:
+	* R/fracdiff.var: finally fix 'FIXME': se.ok
+	* tests/sim*: updates, searching for 'h'  in fracdiff.var(.., h = *)
+
 2022-10-18  Martin Maechler  <[email protected]>
 
 	* DESCRIPTION (Version): 1.5-2:

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: fracdiff
-Version: 1.5-3
-VersionNote: Released 1.5-0 on 2019-12-09, 1.5-1 on 2020-01-20, 1.5-2 on 2022-10-31
-Date: 2024-02-01
+Version: 1.5-4
+VersionNote: Released 1.5-0 on 2019-12-09, 1.5-1 on 2020-01-20, 1.5-2 on 2022-10-31, 1.5-3 on 2024-02-01
+Date: 2024-03-13
 Title: Fractionally Differenced ARIMA aka ARFIMA(P,d,q) Models
 Authors@R: c(person("Martin","Maechler", role=c("aut","cre"), email="[email protected]",
 	   				 comment = c(ORCID = "0000-0002-8685-9910"))

R/fd-methods.R

@@ -19,6 +19,7 @@ logLik.fracdiff <- function(object, ...)
     class(r) <- "logLik"
     r
 }
+## ==> AIC(), BIC() do work; the latter as nobs(logLik(.)) works
 
 print.fracdiff <- function(x, digits = getOption("digits"), ...)
 {

R/fracdiff.R

@@ -1,14 +1,15 @@
+### Patched by Friedrich.Leisch, for use with R, 22.1.1997;  then
+###
+### Copyright 2003--2024 Martin Maechler; fixed, changed enhanced ..
 
 ### Original file:
-### copyright 1991 Department of Statistics, Univeristy of Washington
-
-### Patched by Friedrich.Leisch, for use with R, 22.1.1997
-### fixed & changed by Martin Maechler, since Dec 2003
+### copyright 1991 Department of Statistics, University of Washington
 
 if(getRversion() < "2.15")
 paste0 <- function(...) paste(..., sep="")
 
-.fdcov <- function(x, d, h, nar, nma, hess, fdf.work)
+.fdcov <- function(x, d, h, # <- missing by default
+                   nar, nma, hess, fdf.work)
 {
     npq <- as.integer(nar + nma)
     npq1 <- npq + 1L # integer, too
@@ -150,10 +151,9 @@ fracdiff <- function(x, nar = 0, nma = 0,
     ## Cov == (-H)^{-1} == solve(-H)
 
     ## Note that the following can be "redone" using fracdiff.var() :
-
-    fdc <- .fdcov(x, fdf$d, h,
+    fdc <- .fdcov(x, fdf$d, h, # <- missing by default
                   nar=nar, nma=nma, hess=hess, fdf.work = fdf$w)
-
+    ##==> "vcov" = fdc $ covariance.dpq
     dimnames(hess) <- dimnames(fdc$covariance.dpq)
     hess[1, ] <- fdc$hd
     hess[row(hess) > col(hess)] <- hess[row(hess) < col(hess)]
@@ -172,7 +172,7 @@ fracdiff <- function(x, nar = 0, nma = 0,
 		   n = n,
 		   msg = c(fracdf = fd.msg, fdcov = fdc$msg),
 		   d = fdf$d, ar = fdf$ar, ma = fdf$ma,
-		   covariance.dpq = fdc$covariance.dpq,
+		   covariance.dpq = fdc$covariance.dpq, # == vcov
 		   fnormMin = hstat[2], sigma = sqrt(var.WN),
 		   stderror.dpq	  = if(fdc$se.ok) fdc$stderror.dpq, # else NULL
 		   correlation.dpq= if(fdc$se.ok) fdc$correlation.dpq,
@@ -226,8 +226,7 @@ fracdiff.var <- function(x, fracdiff.out, h)
                           fracdiff.out$ar,
                           rep(0, lwork))),
               info = integer(1))
-## FIXME: should be *automatically* same messages as inside fracdiff() above!
-    fracdiff.out$msg <-
+    msg <-
         if(fdc$info) {
             msg <-
                 switch(fdc$info,
@@ -237,10 +236,14 @@ fracdiff.var <- function(x, fracdiff.out, h)
                        stop("error in gamma function"))
             warning(msg)
         } else "ok"
-    se.ok <- fdc$info != 0 || fdc$info < 3 ## << FIXME -- illogical!!
-    nam <- "d"
-    if(p) nam <- c(nam, paste0("ar", 1:p))
-    if(q) nam <- c(nam, paste0("ma", 1:q))
+    se.ok <- fdc$info %in% 0:2
+
+    if("fdcov" %in% names(fracdiff.out$msg)) # fracdiff(): msg = c(fracdf = fd.msg, fdcov = fdc$msg)
+         fracdiff.out$msg[["fdcov"]] <- msg
+    else fracdiff.out$msg <- msg
+    nam <- c("d",
+             if(p) paste0("ar", 1:p),
+             if(q) paste0("ma", 1:q))
     fracdiff.out$h <- fdc$h
     fracdiff.out$covariance.dpq <- array(fdc$cov, c(npq1,npq1), list(nam,nam))
     fracdiff.out$stderror.dpq    <- if(se.ok) fdc$se # else NULL
@@ -298,3 +301,9 @@ fracdiff.sim <- function(n, ar = NULL, ma = NULL, d, rand.gen = rnorm,
             .Machine$double.eps)[["s"]][ii]
     list(series = y, ar = ar, ma = ma, d = d, mu = mu, n.start = n.start)
 }
+
+## Not exported; used for faster checking, e.g., on CRAN
+doExtras <- function() {
+    interactive() || nzchar(Sys.getenv("R_fracdiff_check_extra")) ||
+        identical("true", unname(Sys.getenv("R_PKG_CHECKING_doExtras")))
+}

man/fd-methods.Rd

@@ -37,13 +37,6 @@
   \item{signif.stars}{logical. If \code{TRUE}, \dQuote{significance stars}
     are printed for each coefficient.}
 }
-% \value{
-%   ~Describe the value returned
-%   If it is a LIST, use
-%   \item{comp1 }{Description of 'comp1'}
-%   \item{comp2 }{Description of 'comp2'}
-%   ...
-% }
 \author{Martin Maechler; Rob Hyndman contributed the
   \code{\link{residuals}()} and \code{\link{fitted}()} methods.}
 \seealso{\code{\link{fracdiff}} to get \code{"fracdiff"} objects,

man/fracdiff.var.Rd

@@ -8,8 +8,9 @@ fracdiff.var(x, fracdiff.out, h)
   \item{x}{a univariate time series or a vector.  Missing values (NAs)
     are not allowed.}
   \item{fracdiff.out}{output from \code{fracdiff} for time series \code{x}.}
-  \item{h}{finite-difference interval for approximating partial
-    derivatives with respect to the \code{d} parameter.}
+  \item{h}{finite-difference interval length (\eqn{ > 0}) for approximating partial
+    derivatives with respect to the \code{d} parameter.  Typically smaller
+    than the one in \code{fracdiff.out}}
 }
 \description{
   Allows the finite-difference interval to be altered for recomputation of the
@@ -27,19 +28,21 @@ fracdiff.var(x, fracdiff.out, h)
 }
 \examples{
 ## Generate a fractionally-differenced ARIMA(1,d,1) model :
-ts.test <- fracdiff.sim(10000, ar = .2, ma = .4, d = .3)
+set.seed(5) # reproducibility; x86_64 Lnx: get warning
+tst <- fracdiff.sim(500, ar = .2, ma = .4, d = .3)$series
 ## estimate the parameters in an ARIMA(1,d,1) model for the simulated series
-fd.out <- fracdiff(ts.test$ser, nar= 1, nma = 1)
+fd.out <- fracdiff(tst, nar= 1, nma = 1) # warning ... maybe change 'h'
+summary(fd.out)## *** Warning ... {has been stored} --> h = 7.512e-6
 
 ## Modify the covariance estimate by changing the finite-difference interval
-(fd.o2 <- fracdiff.var(ts.test$series, fd.out, h = .0001))
+(fd.o2 <- fracdiff.var(tst, fd.out, h = 1e-3))
 ## looks identical as  print(fd.out),
 ## however these (e.g.) differ :
 vcov(fd.out)
 vcov(fd.o2)
 
 ## A case, were the default variance is *clearly* way too small:
-set.seed(1); fdc <- fracdiff(X <- fracdiff.sim(n=100,d=0.25)$series)
+set.seed(1); fdc <- fracdiff(X <- fracdiff.sim(n=100, d=0.25)$series)
 fdc
 # Confidence intervals just based on asymp.normal approx. and std.errors:
 confint(fdc) # ridiculously too narrow

src/fdhess.c

@@ -386,7 +386,7 @@ hesdpq(double *x, double d_, double *hh, double *hd, double *w)
     --w;
 
     /* Function Body */
-    if (*hh <= 0.) {
+    if (*hh <= 0.) { // <==> missing(h) [ = default ] in R's .fdcov()
 	*hh = (fabs(filtfd_.cllf) + 1.) * mauxfd_.epspt5;
     }
     if(*hh > 0.1) *hh = 0.1;

tests/ex.R

@@ -1,7 +1,6 @@
 library(fracdiff)
 
-doExtras <- interactive() # for now
-
+doExtras <- fracdiff:::doExtras()
 .proctime00 <- proc.time()
 
 set.seed(107)
@@ -10,10 +9,11 @@ options(digits = 5)
 ## 1)
 
 x1 <- fracdiff.sim(5000, ar = .2, ma = -.4, d = .3, n.start=0, allow.0 = TRUE)
-(fd1 <- fracdiff(x1$series, nar = 1, nma = 1, dtol = 1e-10))
+summary(fd1 <- fracdiff(x1$series, nar = 1, nma = 1, dtol = 1e-10))
 vcov(fd1)
 logLik(fd1)
-
+stopifnot(all.equal(structure(-7051.5027, df = 4L, nall = 5000L, nobs = 5000L, class = "logLik"),
+                    logLik(fd1)))
 fdCOVcomp <-
     c("h", "covariance.dpq", "stderror.dpq", "correlation.dpq", "hessian.dpq")
 fd1. <- fracdiff.var(x1$series, fd1, h = fd1$h / 2)
@@ -52,7 +52,9 @@ fd1uL <- list(
                            -5875, -5420, 4529,
                             4258,  4529, -5834),
                             3L, 3L, dimnames = dns))
-stopifnot( all.equal(fd1u[fdCOVcomp], fd1uL, tolerance = 2e-4) )
+if(doExtras)
+    print(all.equal(fd1u[fdCOVcomp], fd1uL, tolerance = 0))
+stopifnot(all.equal(fd1u[fdCOVcomp], fd1uL, tolerance = 2e-4) )
 
 ## 2)
 
@@ -107,23 +109,25 @@ sfd2S <- ## dput(sapply(fd2.[fdCOVcomp], signif, digits = 5))
                                  1742.3,  1507.3,-5007.4), 3, dimnames=dns))
 ##
 if(doExtras)
-    print(all.equal(sfd2S, sfd2., tol = 1e-6, countEQ=TRUE)) # 8.7655e-5
+    print(all.equal(sfd2S, sfd2., tol =   0 , countEQ=TRUE)) # 8.7655e-5
 stopifnot(all.equal(sfd2S, sfd2., tol = 2e-4, countEQ=TRUE))
 
 fd2u <- fracdiff.var(x2$series, fd2, h = fd2$h * 8)#-> warning, unable .. corr...
-sd2u <- sapply(fd2u[fdCOVcomp], signif, digits = 4)
+##= no se.ok -->
+fdCOV.0 <- setdiff(fdCOVcomp, c("stderror.dpq", "correlation.dpq"))
+sd2u <- sapply(fd2u[fdCOV.0], signif, digits = 4)
 sd2uS <- list(  ## dput(sapply(sd2u[fdCOVcomp], signif, digits = 5))
     h = 0.0002466,
     covariance.dpq = matrix(c(-0.0003545, 6e-04, 5.724e-05,
                               6e-04, -0.0005003, 5.816e-05,
                               5.724e-05, 5.816e-05, 0.0002371), 3, dimnames=dns),
-    stderror.dpq = c(0, 0, 0.0154),
-    correlation.dpq = matrix(0, 3,3),
+    ## stderror.dpq = c(0, 0, 0.0154),
+    ## correlation.dpq = matrix(0, 3,3),
     hessian.dpq = matrix(c(-3347, -3811, 1742,
                            -3811, -2396, 1507,
                             1742,  1507,-5007), 3, dimnames=dns))
 ##
 if(doExtras)
-    print(all.equal(sd2uS, sd2u, tol = 1e-8, countEQ=TRUE))# 0.000103 (32b Win); T.(64b F30)
+    print(all.equal(sd2uS, sd2u, tol =   0 , countEQ=TRUE))# 0.000103 (32b Win); T.(64b F30)
 stopifnot(all.equal(sd2uS, sd2u, tol = 4e-4, countEQ=TRUE))