improve pcTest() and add tests

GeoBosh · Mar 15, 2020 · a481fd7 · a481fd7
1 parent a4ccfde
commit a481fd7
Show file tree

Hide file tree

Showing 13 changed files with 529 additions and 82 deletions.
diff --git a/.Rbuildignore b/.Rbuildignore
@@ -1,3 +1,4 @@
+.nojekyll
 ^docs/
 .rsync_exclude
 ^_pkgdown.yml
@@ -12,3 +13,4 @@
 ^Experimental
 ^tests/testthat/[.]RData$
 [.]travis[.]yml
+^data-raw$
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -9,7 +9,7 @@ Description: Classes and methods for modelling and simulation of
     <doi:10.1111/j.1467-9892.2009.00617.x>, Boshnakov (1996)
     <doi:10.1111/j.1467-9892.1996.tb00281.x>.
 Version: 0.15-0
-Date: 2020-03-12
+Date: 2020-03-15
 Author: Georgi N. Boshnakov
 Maintainer: Georgi N. Boshnakov <georgi.boshnakov@manchester.ac.uk>
 Depends: R (>= 3.5.0), sarima

diff --git a/NEWS.md b/NEWS.md
@@ -1,5 +1,12 @@
 # Version 0.15.0
 
+- improved `pcTest()`.
+
+- included "nsadata.csv" in the package, currently as internal data for testing.
+
+  Need description and maybe a better name before exporting.  Object `datansa`
+  is the whole dataset, object `nsaauto` is column "AUTOMOTIVEPRODNSA".
+
 - changed slightly some Rd files not rendered well by pkgdown, e.g. move commented
   items out of 'describe' in methods' descriptions. 
 

diff --git a/R/pcTest.R b/R/pcTest.R
@@ -82,10 +82,15 @@ setMethod("pcTest", signature(x = "ANY", nullmodel = "character"),
 
 
 setMethod("pcTest", signature(x = "numeric", nullmodel = "character"),
-          function(x, nullmodel, nseasons, ...){
+          function(x, nullmodel, nseasons, ..., maxlag){
               switch(nullmodel,
                      "wn" = Box.test(x, ...),
                      "piar" = test_piar(x, nseasons, ...),
+                         "pwn" = {
+                              acr <- autocorrelations(x, nseasons = nseasons, maxlag = maxlag)
+                              acrsl <- slMatrix(as.matrix(acr))
+                              pcTest(acrsl, nullmodel, nepoch = floor(length(x)/nseasons), ...)
+                         },
                      # default
                        # do.call(nullmodel, list(x, ...), quote = TRUE)
                        # todo: think about this?
@@ -96,18 +101,30 @@ setMethod("pcTest", signature(x = "numeric", nullmodel = "character"),
 
 
 setMethod("pcTest", signature(x = "PeriodicTimeSeries", nullmodel = "character"),
-          function(x, nullmodel, ...){    # TODO: arg. for a specific column in multivar case?
+          function(x, nullmodel, ..., maxlag){ # TODO: arg. for a specific column in multivar case?
               nseas <- nSeasons(x)
               if(nVariables(x) == 1){
-                      # 2019-04-26 was:  pcTest(coreVector(x), nullmodel, nseas, ...)
-                  pcTest(as(x, "vector"), nullmodel, nseas, ...)
+                  ## 2019-04-26 was:  pcTest(coreVector(x), nullmodel, nseas, ...)
+                  switch(nullmodel,
+                         "pwn" = {
+                              acr <- autocorrelations(x, maxlag)
+                              acrsl <- slMatrix(as.matrix(acr))
+                              pcTest(acrsl, nullmodel, nepoch = nCycles(x), ...)
+                         },
+                         ## default
+                         pcTest(as(x, "vector"), nullmodel, nseas, ...)
+                         )
               }else{ # for now just do the test for each variable separately
                       # 2019-04-26 was:  m <- coreMatrix(x)
-                  m <- as(x, "matrix")
+                      # 2020-03-15 commented out: m <- as(x, "matrix")
                   res <- vector(nVariables(x), mode = "list")
+                  names(res) <- colnames(x)
                   ## TODO: names of variables
                   for(i in seq(along = res))
-                      res[[i]] <- pcTest(m[ , i], nullmodel, nseas, ...)
+                         # 2020-03-15 was: 
+                         #   res[[i]] <- pcTest(m[ , i], nullmodel, nseas, maxlag = maxlag, ...)
+                      res[[i]] <- pcTest(x[[i]], nullmodel, maxlag = maxlag, ...)
+                  res
               }
           }
           )

diff --git a/R/sysdata.rda b/R/sysdata.rda
diff --git a/data-raw/nsadata-csv.R b/data-raw/nsadata-csv.R
@@ -0,0 +1,18 @@
+## code to prepare `nsadata.csv` dataset goes here
+
+## TODO: document and rename, see Lina's thesis
+
+## usethis::use_data("Experimental/nsadata.csv")
+
+datansa <- read.csv("nsadata.csv")
+datansa <- ts(datansa, start = c(1919, 1), frequency = 4)
+
+## now `datansa` is a ts, so this doesn't work:
+## nsaauto_old <- ts(datansa$AUTOMOTIVEPRODNSA[113:328], start = c(1947, 1), frequency = 4)
+##
+## this does: 
+##    tmp <- datansa[ , "AUTOMOTIVEPRODNSA"][113:328]
+##    nsaauto_old <- ts(tmp, start = c(1947, 1), frequency = 4)
+## but it can be done also this way:
+nsaauto <- window(datansa[ , "AUTOMOTIVEPRODNSA"], start = c(1947, 1))
+## identical(nsaauto, nsaauto_old) # TRUE
diff --git a/data-raw/nsadata.csv b/data-raw/nsadata.csv
diff --git a/docs/news/index.html b/docs/news/index.html
diff --git a/inst/REFERENCES.bib b/inst/REFERENCES.bib
@@ -1,3 +1,14 @@
+@article{francq2011asymptotic,
+  title={Asymptotic properties of weighted least squares estimation in weak parma models},
+  author={Francq, C. and Roy, R. and Saidi, A.},
+  journal={Journal of Time Series Analysis},
+  volume={32},
+  number={6},
+  pages={699--723},
+  year={2011},
+  publisher={Wiley Online Library}
+}
+
 
 @article{boswijk1996unit,
   title = {Unit roots in periodic autoregressions},

diff --git a/man/mC.ss.Rd b/man/mC.ss.Rd
@@ -158,75 +158,3 @@ xxcoz4a <- mC.ss(spec.coz4)
 % R documentation directory.
 \keyword{ pcts }
 \keyword{ torevise }
-
-% > xxco.1 <- mC.ss(m1.new, generators = TRUE)
-%
-% > datansa <- read.csv("nsadata.csv")
-% > nsaauto <- ts(datansa$AUTOMOTIVEPRODNSA[113:328], start=c(1947, 1), frequency=4)
-%
-% > res.xxco.1 <- xxco.1$env$minimBB(nsaauto, control=list(maxit=1000))
-%
-% condlik is:  32.85753	persd is:  16.96771 10.40725 3.567698 7.426556
-% iter:  0  f-value:  32.85753  pgrad:  14.83674
-% iter:  10  f-value:  30.21297  pgrad:  0.0007615952
-%   Successful convergence.
-%
-% > res.xxco.1$value
-% [1] 30.21297
-% > res.xxco.1$par
-%      co.r1      co.r2      co.r3      co.r4
-% -0.4069477 -0.5093360 -0.6026860 -0.5174826
-% > res.xxco.1
-% $par
-%      co.r1      co.r2      co.r3      co.r4
-% -0.4069477 -0.5093360 -0.6026860 -0.5174826
-%
-% $value
-% [1] 30.21297
-%
-% $gradient
-% [1] 9.023893e-06
-%
-% $fn.reduction
-% [1] 2.644559
-%
-% $iter
-% [1] 14
-%
-% $feval
-% [1] 16
-%
-% $convergence
-% [1] 0
-%
-% $message
-% [1] "Successful convergence"
-%
-% $cpar
-% method      M
-%      2     50
-%
-% > with(xxco.1$env, model)
-% $period
-% [1] 4
-%
-% $p
-% [1] 5
-%
-% $q
-% [1] 0
-%
-% $phi
-%           [,1]          [,2]          [,3] [,4]       [,5]
-% [1,] 1.1646497 -1.165471e-16 -4.254923e-17    1 -1.1646497
-% [2,] 0.8451102 -2.220446e-16 -5.456035e-17    1 -0.8451102
-% [3,] 0.7989768  0.000000e+00  2.220446e-16    1 -0.7989768
-% [4,] 1.2716195 -1.110223e-16 -6.058867e-17    1 -1.2716195
-%
-% > with(xxco.1$env, zapsmall(model$phi))
-%           [,1] [,2] [,3] [,4]       [,5]
-% [1,] 1.1646497    0    0    1 -1.1646497
-% [2,] 0.8451102    0    0    1 -0.8451102
-% [3,] 0.7989768    0    0    1 -0.7989768
-% [4,] 1.2716195    0    0    1 -1.2716195
-
diff --git a/tests/testthat/test-fitPM.R b/tests/testthat/test-fitPM.R
@@ -73,3 +73,113 @@ test_that("test fitPM()",
     fitPM(pipfm, perunit)
 
 })
+
+
+test_that("test mC.ss() works",
+{
+    ## examples from mC.ss.Rd
+# test0 roots
+spec.coz2 <- mcompanion::mcSpec(dim = 5, mo = 4, root1 = c(1,1), order = rep(2,4))
+spec.coz2
+xxcoz2a <- mC.ss(spec.coz2)
+
+## test0 roots
+spec.coz4 <- mcompanion::mcSpec(dim = 5, mo = 4, root1 = c(1,1), order = rep(3,4))
+xxcoz4a <- mC.ss(spec.coz4)
+
+    ## excerpt from 
+    ## ~/Documents/Rwork/pctsExperiments/Rsessions/combined upto 2013-12-31 17h36m.Rhistory
+spec.co2 <- mcompanion::mcSpec(dim = 5, mo = 4, siorder = 1)
+tmp2 <- mC.ss(spec.co2)
+## only two iters for testthat
+expect_output(mc.res1ssenv2b <- tmp2$env$minimBB(nsaauto, control=list(maxit=2)))
+
+expect_output(tmp2$env$minimBB(nsaauto, control=list(maxit=2)))
+expect_output(tmp2$env$minimBBlu(nsaauto, control=list(maxit=2)))
+expect_output(tmp2$env$minimBB(nsaauto, control=list(maxit=2), CONDLIK = FALSE))
+tmp2$env$minim(nsaauto, control=list(maxit=2))
+tmp2$env$minim(nsaauto, control=list(maxit=2), CONDLIK = FALSE)
+expect_output(tmp2$env$minimBB(nsaauto, control=list(maxit=2), CONDLIK = FALSE))
+mC.ss(spec.co2, generators = TRUE)
+
+    tmp2$env$mcparam2optparam()
+    tmp2$env$mcsigma2(nsaauto)
+    tmp2$env$mcsigma2(nsaauto, tmp2$env$mcparam2optparam())
+
+mC.ss(spec.co2, init = tmp2$env$mcparam2optparam())    
+
+    ## this chunk was commented out in mC.ss.Rd, old testing with it.
+##  > xxco.1 <- mC.ss(m1.new, generators = TRUE)
+## 
+##  > datansa <- read.csv("nsadata.csv")
+##  > nsaauto <- ts(datansa$AUTOMOTIVEPRODNSA[113:328], start=c(1947, 1), frequency=4)
+## 
+##  > res.xxco.1 <- xxco.1$env$minimBB(nsaauto, control=list(maxit=1000))
+## 
+##  condlik is:  32.85753	persd is:  16.96771 10.40725 3.567698 7.426556
+##  iter:  0  f-value:  32.85753  pgrad:  14.83674
+##  iter:  10  f-value:  30.21297  pgrad:  0.0007615952
+##    Successful convergence.
+## 
+##  > res.xxco.1$value
+##  [1] 30.21297
+##  > res.xxco.1$par
+##       co.r1      co.r2      co.r3      co.r4
+##  -0.4069477 -0.5093360 -0.6026860 -0.5174826
+##  > res.xxco.1
+##  $par
+##       co.r1      co.r2      co.r3      co.r4
+##  -0.4069477 -0.5093360 -0.6026860 -0.5174826
+## 
+##  $value
+##  [1] 30.21297
+## 
+##  $gradient
+##  [1] 9.023893e-06
+## 
+##  $fn.reduction
+##  [1] 2.644559
+## 
+##  $iter
+##  [1] 14
+## 
+##  $feval
+##  [1] 16
+## 
+##  $convergence
+##  [1] 0
+## 
+##  $message
+##  [1] "Successful convergence"
+## 
+##  $cpar
+##  method      M
+##       2     50
+## 
+##  > with(xxco.1$env, model)
+##  $period
+##  [1] 4
+## 
+##  $p
+##  [1] 5
+## 
+##  $q
+##  [1] 0
+## 
+##  $phi
+##            [,1]          [,2]          [,3] [,4]       [,5]
+##  [1,] 1.1646497 -1.165471e-16 -4.254923e-17    1 -1.1646497
+##  [2,] 0.8451102 -2.220446e-16 -5.456035e-17    1 -0.8451102
+##  [3,] 0.7989768  0.000000e+00  2.220446e-16    1 -0.7989768
+##  [4,] 1.2716195 -1.110223e-16 -6.058867e-17    1 -1.2716195
+## 
+##  > with(xxco.1$env, zapsmall(model$phi))
+##            [,1] [,2] [,3] [,4]       [,5]
+##  [1,] 1.1646497    0    0    1 -1.1646497
+##  [2,] 0.8451102    0    0    1 -0.8451102
+##  [3,] 0.7989768    0    0    1 -0.7989768
+##  [4,] 1.2716195    0    0    1 -1.2716195
+
+
+
+})
diff --git a/tests/testthat/test-pcarma.R b/tests/testthat/test-pcarma.R
@@ -13,9 +13,9 @@ pcarma_unvec(list(p = 2, q = 0, period = 2, param = coef3))
 ## actually, the model is PAR(1,2):
 s3a <- pcarma_param_system(pc3, NULL, NULL, c(1, 2), 0, 2)
 coef3a <- solve(s3a$A, s3a$b)
-pcarma_unvec(list(p = c(1,2), q = 0, period = 2, param = coef3a))
-
-
+coef3a_more <- pcarma_unvec(list(p = c(1,2), q = 0, period = 2, param = coef3a))
+coef3a_vec <- pcarma_tovec(coef3a_more)
+    
 ## prepare test parameters for a PAR(2) model with period=2.
 ##   (rounded to 6 digits from the above example.
 m1 <- rbind(c(1, 0.81, 0), c(1, 0.4972376, 0.4972376) )

diff --git a/tests/testthat/test-pclspiar.R b/tests/testthat/test-pclspiar.R
@@ -36,4 +36,24 @@ test_that("pclspiar() is ok", {
     ##     LRurpar.test(cu, list(regular = c(0,0,0), seasonal = c(1,0), regvar = 0), p = 1)
     ## }
 
+
+    pcTest(pcts(nsaauto), "wn")
+    pcTest(pcts(nsaauto), "piar", p = 1)
+
+    ## tmpslMat <- slMatrix(rnorm(32), period = 4, maxlag = 7)
+
+    acr <- autocorrelations(pcts(nsaauto), maxlag = 7)
+    acrsl <- slMatrix(as.matrix(acr))
+    pcTest(acrsl, "pwn", nepoch = nCycles(pcts(nsaauto)))
+    pcTest(acrsl, "periodicity", nepoch = nCycles(pcts(nsaauto)))
+
+    pcTest(pcts(nsaauto), "pwn", maxlag = 4)
+    pcTest(as.numeric(nsaauto), "pwn", maxlag = 4, nseasons = 4)
+    pcTest(acrsl, "periodicity", nepoch = nCycles(pcts(nsaauto)))
+    pcTest(pcts(nsaauto), "wn")
+    pcTest(pcts(nsaauto), "wn", lag = 3)
+    ## pcTest(matrix(nsaauto, nrow = 4), "wn")
+
+    pcTest(pcts(datansa), "pwn", maxlag = 4)
+
 })