context("fitPM") test_that("test fitPM()", { set.seed(1234) x <- arima.sim(list(ar = 0.9), n = 1000) mx <- matrix(x, nrow = 4) x_pcts <- pcts(as.numeric(x), nseasons = 4) expect_error(fitPM(c(1.5, 2, 3, 1), x), "The PAR orders must be non-negative integer numbers") expect_error(fitPM("dummy", x), "unable to find an inherited method for function 'fitPM'") expect_error(fitPM(c(3,2,2,2), mx), "multivariate PAR fitting not implemented yet") proba1 <- fitPM(c(3, 2, 2, 2), as.numeric(mx)) ## in package Matrix v1.5.0 slot asyCov is "dsTMatrix" (i.e. sparse symmetric), while in ## previous version of Matrix it was "dgTMatrix" (sparse but not symmetric) causing check ## below to fail with 1.5.0. ## ## I recomputed the value saved in "proba1.RDS" with v1.5.0 so that it works with it and ## for now skip the check with older versions of Matrix. ## ## With Matrix v1.5.2 the value of asyCov returned to the pre 1.5.0 class, due to ## multiplication of 'matrix' by a Permutation Matrix now giving Matrix instead of ## 'matrix', which changes the class of some arguments to .cbind(). Need to do the ## comparison below class independent. See the email exchange with Mikael Jagan on ## 2022-10-19. ## ## Commenting this test for now. ## if(packageVersion("Matrix") >= '1.5.0') ## expect_equal_to_reference(proba1, "proba1.RDS") expect_output(show(proba1)) expect_output(summary(proba1)) expect_error(fitPM(2, mx), ## "unable to find an inherited method for function [.]nSeasons[.] for signature [.]\"matrix\"[.]" "unable to find an inherited method for function" ) expect_identical(fitPM(2, x_pcts), fitPM(c(2, 2, 2, 2), as.numeric(mx))) eps.proba1 <- residuals(proba1) expect_identical(eps.proba1, pcts:::.whiten(proba1)) eps.proba1[1:2] <- NA dim(eps.proba1) <- dim(mx) pc.cconesidedsum(mx, eps.proba1, maxlag = 4) ## estimate h_{t,i}, see Boshnakov (1996) pc.hat.h(mx, eps.proba1, maxlag = 4) ## replacing with equivalent code using the new Fraser2017 ## data(Fraser, package = "pear") Fraser <- window(Fraser2017, start = c(1912, 3), end = c(1990, 12)) logFraser <- log(Fraser) ## TODO: for now I need whole years; ## !!! However note the typo 'logfraser', the following use 'logFraser'! logfraser <- ts(logFraser[1:936], frequency = 12) #### co1_pear <- pear::pear(logFraser, 1)[["phi"]] # fitPM(as.numeric(logFraser), order = rep(1, 12), period = 12, seasonof1st = 3) az1 <- fitPM(model = rep(1, 12), as.numeric(logFraser), seasonof1st = 3) az2 <- fitPM(model = rep(1, 12), as.numeric(logFraser)) #### expect_true(all.equal(as.vector(az1@ar@coef[ , 1]), as.vector(co1_pear[ , 1]))) #### expect_true(all.equal(as.vector(az2@ar@coef[ , 1]), as.vector(co1_pear[ , 1])[c(3:12, 1:2)])) ## pcfr2 <- pcts(dataFranses1996[ , 2 ]) pcfr23 <- pcts(dataFranses1996[ , 2:3]) expect_error(fitPM(model = rep(1, 4), pcfr23), "Multivariate case not implemented yet") ## fitPM(model = rep(1, 4), pcfr23[1]) # tests the method for PeriodicMTS ([] keep MTS class) ## fitPM(model = rep(1, 4), pcfr23[[1]]) # tests the method for PeriodicTS ('[[' drops the 'M') expect_identical(fitPM(model = rep(1, 4), pcfr23[1]), fitPM(model = rep(1, 4), pcfr23[[1]])) x <- arima.sim(list(ar = 0.9), n = 960) pcx <- pcts(x, nseasons = 4) mx <- matrix(x, nrow = 4) pfm1 <- PeriodicArModel(matrix(1:12, nrow = 4), order = rep(3,4), sigma2 = 1) sipfm1 <- new("SiPeriodicArModel", iorder = 1, siorder = 1, pcmodel = pfm1) fitPM(sipfm1, mx) expect_output(show(sipfm1)) d4piar2 <- rbind(c(1,0.5,-0.06), c(1, 0.6, -0.08), c(1, 0.7, -0.1), c(1, 0.2, 0.15)) picoef1 <- c(0.8, 1.25, 2, 0.5) parcoef1 <- d4piar2[, 2:3] coef1 <- pi1ar2par(picoef1, parcoef1) tmpval <- PeriodicArModel(parcoef1) ##pipfm <- PiParModel(piorder = 1, picoef = picoef1, par = tmpval) pipfm <- new("PiPeriodicArModel", piorder = 1, picoef = matrix(picoef1, ncol = 1), pcmodel = tmpval) expect_output(show(pipfm)) perunit <- sim_pc(list(phi = coef1, p = 3, q = 0, period = 4),500) fitPM(pipfm, perunit) ## temporary proba1x <- new("FittedPeriodicArmaModel", as(proba1, "PeriodicArmaModel"), theTS = proba1@theTS, ns = proba1@ns, asyCov = proba1@asyCov) expect_identical(residuals(proba1x), residuals(proba1)) expect_error(as_pcarma_list(1:10), "unable to find an inherited method for function 'as_pcarma_list'") expect_output(show(proba1x)) fitted(proba1x) predict(proba1x, 1) predict(proba1x, 8) n <- 100 x <- arima.sim(list(ar=0.9), n = n) proba1 <- fitPM(c(3,2,2,2), x) meancovmat(proba1, n/10) meancovmat(proba1, n/10, cor = TRUE) meancovmat(proba1, n/10, result = "") meancovmat(proba1, n/10, cor = TRUE, result = "") meanvarcheck(proba1, 100) }) test_that("test mC.ss() works", { pcts_exdata() ## 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 set.seed(1234) ## This prints something like: ## condlik is: 18.09375 persd is: 1.651785 1.714789 3.041003 1.577415 fitPM(spec.coz4, rnorm(100), control = list(maxit = 1)) })