Hyndman's residuals() and fitted(); using Authors@R ; more examples

Author: mmaechler

.Rbuildignore

@@ -6,7 +6,9 @@ Done
 .*\.tgz
 .*\.mail
 .*-manual
+qed.*\.pdf
 codes_for_.*
+^filters\.R
 simul-compare.*
 tests/.*linux
 tests/windows.*

ChangeLog

@@ -1,3 +1,15 @@
+2019-12-09  Martin Maechler  <[email protected]>
+
+	* DESCRIPTION (Authors@R): using new format; ready for CRAN
+
+	* R/fd-methods.R (summary.fracdiff): renamed '$ coefficients' such
+	that 'coef()' works.
+
+2019-12-09  Rob Hyndman  <[email protected]>
+
+	* R/fracdiff.R, fd-methods.R (fracdiff): provide 'residuals()' and
+	'fitted'()' methods.
+
 2018-09-10  Martin Maechler  <[email protected]>
 
 	* R/diffseries.R (diffseries): Finally found bug in the fft() based

DESCRIPTION

@@ -1,15 +1,25 @@
 Package: fracdiff
 Version: 1.5-0
 VersionNote: Released 1.4-2 on 2012-12-01 on CRAN
-Date: 2018-09-10
+Date: 2019-12-09
 Title: Fractionally Differenced ARIMA aka ARFIMA(P,d,q) Models
-Author: S original by Chris Fraley, U.Washington, Seattle.
-   R port by Fritz Leisch at TU Wien; since 2003-12: Martin Maechler;
-   fdGPH(), fdSperio(), etc by Valderio Reisen and Artur Lemonte.
-Maintainer: Martin Maechler <[email protected]>
+Authors@R: c(person("Martin","Maechler", role=c("aut","cre"), email="[email protected]",
+	   				 comment = c(ORCID = "0000-0002-8685-9910"))
+ , person("Chris", "Fraley", role=c("ctb","cph"), comment = "S original; Fortran code")
+ , person("Friedrich", "Leisch", role = "ctb",
+                            comment = c("R port", ORCID = "0000-0001-7278-1983"))
+ , person("Valderio", "Reisen", role="ctb", comment = "fdGPH() & fdSperio()")
+ , person("Artur",  "Lemonte",  role="ctb", comment = "fdGPH() & fdSperio()")
+ , person("Rob", "Hyndman", email="[email protected]", role="ctb",
+     	                    comment = c("residuals() & fitted()", ORCID = "0000-0002-2140-5352"))
+ )
 Description: Maximum likelihood estimation of the parameters of a fractionally
-   differenced ARIMA(p,d,q) model (Haslett and Raftery, Appl.Statistics, 1989).
+   differenced ARIMA(p,d,q) model (Haslett and Raftery, Appl.Statistics, 1989);
+   including inference and basic methods.  Some alternative algorithms to estimate "H".
 Imports: stats
-Suggests: longmemo, urca
+Suggests: longmemo, forecast, urca
 License: GPL (>= 2)
+URL: https://github.com/mmaechler/fracdiff
+BugReports: https://github.com/mmaechler/fracdiff/issues
 Encoding: UTF-8
+NeedsCompilation: yes

NAMESPACE

@@ -1,14 +1,16 @@
 useDynLib(fracdiff, .registration=TRUE)
 
 importFrom("stats",
-           AIC, as.ts
+           AIC, arima, as.ts
          , fft, nextn
          , lm.fit
          , coef, logLik, resid
          , na.fail
          , printCoefmat
          , pnorm, qnorm, rnorm
-         , symnum)
+         , symnum
+         , tsp, "tsp<-"
+           )
 
 export("fracdiff", "fracdiff.sim", "fracdiff.var",
        "diffseries",

R/fd-methods.R

@@ -5,8 +5,11 @@ coef.fracdiff <- function(object, ...) unlist(object[c("d", "ar", "ma")])
 
 vcov.fracdiff      <- function(object, ...) object$covariance.dpq
 
-residuals.fracdiff <- function(object, ...) .NotYetImplemented()
-fitted.fracdiff    <- function(object, ...) .NotYetImplemented()
+## Lines added by RJH. 9 Dec 2019
+residuals.fracdiff <- function(object, ...) object$residuals
+
+fitted.fracdiff    <- function(object, ...) object$fitted
+## End of RJH addition
 
 logLik.fracdiff <- function(object, ...)
 {
@@ -40,7 +43,7 @@ summary.fracdiff <- function(object, symbolic.cor = FALSE, ...)
     cf <- cbind("Estimate" = cf,
                 "Std. Error"= se, "z value" = cf / se,
                 "Pr(>|z|)" = 2 * pnorm(-abs(cf / se)))
-    object$coef <- cf
+    object$coefficients <- cf # 'long name' such that coef(.) works
     logl <- logLik(object)
     object$df <- attr(logl, "df")
     object$aic <- AIC(logl)

R/fracdiff.R

@@ -94,6 +94,7 @@ fracdiff <- function(x, nar = 0, nma = 0,
     if(is.null(dtol))
         dtol <- .Machine$double.eps^0.25 # ~ 1.22e-4
     ## if dtol < 0: the fortran code will choose defaults
+    tspx <- tsp(x) # Added by RJH. 9 Dec 2019
     x <- as.double(x)
 
     ## this also initializes "common blocks" that are used in .C(.) calls :
@@ -159,16 +160,26 @@ fracdiff <- function(x, nar = 0, nma = 0,
 
     hstat <- fdf[["hood.etc"]]
     var.WN <- hstat[3]
+
+    ## Following lines added by RJH. 9 Dec 2019
+    diffx <- diffseries(x, d = fdf$d)
+    armafit <- arima(diffx, order = c(length(fdf$ar), 0L, length(fdf$ma)),
+                     include.mean = FALSE, fixed = c(fdf$ar, -fdf$ma))
+    res <- armafit$residuals
+    tsp(res) <- tspx
+
     structure(list(log.likelihood = hstat[1],
-                   n = n,
+		   n = n,
 		   msg = c(fracdf = fd.msg, fdcov = fdc$msg),
 		   d = fdf$d, ar = fdf$ar, ma = fdf$ma,
 		   covariance.dpq = fdc$covariance.dpq,
 		   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,
 		   h = fdc$h, d.tol = fdf$dtol, M = M, hessian.dpq = hess,
-		   length.w = lenw, call = cl),
+		   length.w = lenw,
+		   residuals = res, fitted = x - res, ## by RJH
+		   call = cl),
 	      class = "fracdiff")
 }
 

man/fd-methods.Rd

@@ -22,8 +22,8 @@
 \method{print}{summary.fracdiff}(x, digits = max(3, getOption("digits") - 3),
         correlation = FALSE, symbolic.cor = x$symbolic.cor,
         signif.stars = getOption("show.signif.stars"), \dots)
-% \method{fitted}{fracdiff}(object, \dots)
-% \method{residuals}{fracdiff}(object, \dots)
+\method{fitted}{fracdiff}(object, \dots)
+\method{residuals}{fracdiff}(object, \dots)
 \method{vcov}{fracdiff}(object, \dots)
 }
 \arguments{
@@ -37,17 +37,15 @@
   \item{signif.stars}{logical. If \code{TRUE}, \dQuote{significance stars}
     are printed for each coefficient.}
 }
-% \details{
-%   ~~ If necessary, more details than the description above ~~
-% }
 % \value{
 %   ~Describe the value returned
 %   If it is a LIST, use
 %   \item{comp1 }{Description of 'comp1'}
 %   \item{comp2 }{Description of 'comp2'}
 %   ...
 % }
-\author{Martin Maechler}
+\author{Martin Maechler; Rob Hyndman contributed the
+  \code{\link{residuals}()} and \code{\link{fitted}()} methods.}
 \seealso{\code{\link{fracdiff}} to get \code{"fracdiff"} objects,
   \code{\link{confint.fracdiff}} for the \code{\link{confint}} method;
   further, \code{\link{fracdiff.var}}.
@@ -63,6 +61,9 @@ smFD <- summary(modFD)
 smFD
 coef(smFD) # gives the whole table
 AIC(modFD) # AIC works because of the logLik() method
+stopifnot(exprs = {
+
+})
 }
 \keyword{print}
 \keyword{models}

tests/ex.Rout-64b

@@ -1,8 +1,7 @@
 
-R version 2.13.1 Patched (2011-08-08 r56671)
-Copyright (C) 2011 The R Foundation for Statistical Computing
-ISBN 3-900051-07-0
-Platform: x86_64-unknown-linux-gnu (64-bit)
+R version 3.6.2 RC (2019-12-06 r77534) -- "Dark and Stormy Night"
+Copyright (C) 2019 The R Foundation for Statistical Computing
+Platform: x86_64-pc-linux-gnu (64-bit)
 
 R is free software and comes with ABSOLUTELY NO WARRANTY.
 You are welcome to redistribute it under certain conditions.
@@ -40,7 +39,7 @@ a list with components:
  [5] "ar"              "ma"              "covariance.dpq"  "fnormMin"       
  [9] "sigma"           "stderror.dpq"    "correlation.dpq" "h"              
 [13] "d.tol"           "M"               "hessian.dpq"     "length.w"       
-[17] "call"           
+[17] "residuals"       "fitted"          "call"           
 > vcov(fd1)
               d         ar1         ma1
 d    0.00059661 -0.00080522 -0.00018971
@@ -120,7 +119,7 @@ d    0.3374     0.0205    16.5   <2e-16 ***
 ar   0.7710     0.0414    18.6   <2e-16 ***
 ma  -0.3810     0.0160   -23.8   <2e-16 ***
 ---
-Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
+Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 sigma[eps] = 1.0088 
 [d.tol = 0.0001221, M = 100, h = 3.082e-05]
 Log likelihood: -2.92e+03 ==> AIC = 5856.5 [4 deg.freedom]
@@ -183,5 +182,8 @@ ma1  1742  1507 -5007
 > 
 > ## Last Line:
 > cat('Time elapsed: ', proc.time() - .proctime00,'\n')
-Time elapsed:  0.21 0.003 0.213 0 0 
+Time elapsed:  0.065 0.002 0.067 0 0 
 > 
+> proc.time()
+   user  system elapsed 
+  0.164   0.021   0.225