tests/testthat/test-Cpp_arma_related_functions.R

In deseats: Data-Driven Locally Weighted Regression for Trend and Seasonality in TS

Xt <- c(12.8433856964746, 11.6161505776814, 11.3094186367253, 10.7116011400713, 10.1862273446052, 10.0708742587288, 8.98039903644311, 10.2107197146335, 11.6497704403272, 11.0233889301354, 12.3311187036587, 13.6151547745729, 13.3787689699471, 14.0427281985119, 15.0856060570139, 15.7505720653137, 16.0548312406384, 16.0666196490345, 15.3916808705163, 14.3672580903331, 13.4536283303508, 12.4625973518189, 12.0985473787934, 11.0296281745773, 13.1869975117585, 15.4533382392052, 14.2282091833928, 12.9253073014742, 12.3096173397592, 13.0933695793561, 13.4019721982456, 13.5230145195308, 13.4942226653202, 13.3188120307585, 14.6003354768997, 14.6087649901778, 15.5597206047244, 14.0312869212916, 13.5621381924251, 13.8396578533602, 13.8956291880184, 14.1445506987448, 13.518681777392, 12.8087081339503, 11.7109166188234, 10.5465627151473, 11.0802269287696, 12.2861329494702, 12.8074021806744, 13.8439511077171, 16.131984194969, 15.7564001132353, 12.4047956205956, 12.4645110210588, 12.4488262143899, 11.6015693523144, 12.7326656080831, 13.0708786034555, 11.7165108066388, 11.702502188269, 12.0735119092239, 12.3177768411371, 12.9406062828722, 12.8104045667813, 13.3563297293695, 13.3811325229993, 13.5183314491059, 14.6555620094813, 15.1724944040454, 14.4424052427174, 15.0255254955877, 16.0157755120991, 16.1861917459585, 15.7732644370364, 14.4627448974418, 15.1006208387196, 14.5805012242216, 16.0235540563981, 17.7727649621722, 16.9015404275126, 14.6168072194018, 12.9381480397286, 12.8883840551169, 12.8336188146977, 12.4426446316595, 11.5144613751978, 11.469285954197, 11.1510818360764, 9.67536311393678, 9.48855669117392, 11.3008730327386, 11.724762464839, 12.4576245663082, 11.3364366492391, 11.0010023347982, 12.1223278555399, 12.9764864630164, 13.2697933993438, 12.6309874683097, 11.6016338776306, 10.6021627507197, 10.9570961273014, 10.7113159734763, 10.4636264306403, 10.7531675873996, 13.1730332203128, 13.4211043551396, 13.2341846700349, 13.3895834047305, 12.3798321739631, 12.0436461185855, 13.7694490983714, 14.7423938511509, 14.5166606443382, 13.6897820796936, 11.2483518180752, 11.9085598752965, 11.3710204484898, 11.8993417893578, 14.5386165772066, 13.6176902276816, 13.4422209142659, 13.428249318271, 11.626531618943, 9.7678323999056, 8.40861267946845, 8.48435426667032, 8.14183914058817, 9.5863263788896, 13.0844957922363, 12.8521710487931, 13.0707433436127, 14.2424320052674, 14.6164261031378, 13.3192229533945, 12.629815156607, 12.4307489711393, 13.0543161454329, 12.9770773123094, 13.8047081272503, 13.707201524651, 14.3904621039515, 13.4805021988897, 11.5816077911316, 14.6480877663555, 15.3873238157596, 14.7909254684922, 15.0693312977488, 14.3130361197323, 14.2240970152719, 14.5142634602322, 14.0590861782134, 13.7289631044538, 13.5388534361419, 15.5186156774463, 15.1159208274476, 13.0444802640023, 12.5509264504353, 13.0450234643334, 13.6532802786651, 13.2471926018096, 11.91052602094, 12.9988903237619, 13.3111654267245, 12.207477742077, 11.8061984191861, 11.8790232146112, 13.2566485927672, 13.849351451991, 14.4209684476769, 13.9024127887938, 13.6264882531511, 13.2510785205058, 13.1379955367822, 12.2785252818568, 10.766484484983, 12.8083517575721, 14.4951617463469, 13.1256962726541, 11.881185793365, 10.7557670124179, 13.0651072152776, 15.4873932015697, 15.1622366726192, 14.9541872017174, 14.2551881228993, 13.2252148174545, 12.1609240870601, 11.4537281745575, 13.3012085557972, 14.0122715976409, 13.8308212922335, 13.9148849523763, 15.0247390599483, 14.5474914502617, 12.8627049493188, 14.1133368217303, 14.2060210245075, 12.9305671387628, 11.3835908932807, 9.99566134353704, 9.71063160001905, 11.0584411975309, 13.0897498857846, 14.3553373433497, 13.9688730201673, 13.5861324029541, 12.7655711870801, 11.796383437542, 12.7135098244693, 12.2561895480134, 13.9568303598001, 14.5871086406405, 14.3103550780946, 13.3159902632413, 12.3075406421358, 10.9410679723082, 10.7473141260904, 11.7806395990014, 12.7102432659246, 12.2077859711731, 11.2960599654398, 10.940732272038, 12.8334739011015, 12.5241584450285, 11.9492075068069, 14.1072812792536, 14.6553029209882, 12.7783838982059, 11.5301789415594, 12.1542005203594, 12.3050674820343, 11.7938694380393, 11.5608531907242, 11.9411352635819, 13.9325517792228, 14.39347551889, 15.0419498033051, 15.6305187260399, 15.073100270547, 12.8914482387393, 11.6664454481197, 11.346082701769, 11.6744770974306, 13.425909862249, 16.2713152162388, 18.0140411972665, 17.1804684032495, 14.1332250045525, 12.5562139413362, 12.6141670499728, 13.649533550585, 14.8697079559275, 15.5043689828079, 13.7299386528743, 13.6929629277456, 13.4933383558114, 13.3268609410156, 13.4207551859653, 13.7899709837331, 13.8217380539245, 11.9500790315746, 12.1978137380828, 13.2227441673944, 13.1132825760727, 13.0946165016584, 13.2893424520527, 13.5205785076821, 15.1457608556355, 15.1334053030001, 14.581923906518, 15.389763299574, 16.3412072696133, 17.0411987107981, 15.8540026516893, 16.7415505859456, 16.7170276462385, 17.5321216325389, 15.7850083790999, 14.3250381087088, 15.2179330402425, 14.4922682264316, 12.964277696673, 11.6901111999553, 10.5912952088265, 10.3804904733288, 11.3089427376989, 9.99763900494803, 10.1332224878756, 12.3659210046156, 15.2420562446796, 16.9194023725966, 17.1809444207138, 14.6083180772833, 12.6909507465344, 12.1608657119881, 12.9701066728572, 13.2551652883747, 11.8971314068003, 13.3131033788242)
n <- length(Xt)
ln <- log(n)

p <- 1
q <- 0
fit1_a <- armaCpp(Xt, p = p, q = q, armamean = TRUE)
fit1_b <- stats::arima(Xt, order = c(p, 0, q))
BIC1_a <- BICarmaCpp(Xt, p = p, q = q, armamean = TRUE)
BIC1_b <- -2 * fit1_b$loglik + ln * (p + q)

p <- 0
q <- 1
fit2_a <- armaCpp(Xt, p = p, q = q, armamean = TRUE)
fit2_b <- stats::arima(Xt, order = c(p, 0, q))
BIC2_a <- BICarmaCpp(Xt, p = p, q = q, armamean = TRUE)
BIC2_b <- -2 * fit2_b$loglik + ln * (p + q)

p <- 0
q <- 0
fit3_a <- armaCpp(Xt, p = p, q = q, armamean = TRUE)
fit3_b <- stats::arima(Xt, order = c(p, 0, q))
BIC3_a <- BICarmaCpp(Xt, p = p, q = q, armamean = TRUE)
BIC3_b <- -2 * fit3_b$loglik + ln * (p + q)

p <- 2
q <- 0
fit4_a <- armaCpp(Xt, p = p, q = q, armamean = TRUE)
fit4_b <- stats::arima(Xt, order = c(p, 0, q))
BIC4_a <- BICarmaCpp(Xt, p = p, q = q, armamean = TRUE)
BIC4_b <- -2 * fit4_b$loglik + ln * (p + q)

p <- 0
q <- 2
fit5_a <- armaCpp(Xt, p = p, q = q, armamean = TRUE)
fit5_b <- stats::arima(Xt, order = c(p, 0, q))
BIC5_a <- BICarmaCpp(Xt, p = p, q = q, armamean = TRUE)
BIC5_b <- -2 * fit5_b$loglik + ln * (p + q)

p <- 2
q <- 2
fit6_a <- armaCpp(Xt, p = p, q = q, armamean = TRUE)
fit6_b <- stats::arima(Xt, order = c(p, 0, q))
BIC6_a <- BICarmaCpp(Xt, p = p, q = q, armamean = TRUE)
BIC6_b <- -2 * fit6_b$loglik + ln * (p + q)

### Basic ARMA fitting with different orders using 'armaCpp'

test_that("armaCpp works as intended", {
  
  expect_equal(fit1_a$coef, fit1_b$coef)
  expect_equal(fit1_a$residuals, fit1_b$residuals)
  
  expect_equal(fit2_a$coef, fit2_b$coef)
  expect_equal(fit2_a$residuals, fit2_b$residuals)
  
  expect_equal(fit3_a$coef, fit3_b$coef)
  expect_equal(fit3_a$residuals, fit3_b$residuals)  
  
  expect_equal(fit4_a$coef, fit4_b$coef)
  expect_equal(fit4_a$residuals, fit4_b$residuals)
  
  expect_equal(fit5_a$coef, fit5_b$coef)
  expect_equal(fit5_a$residuals, fit5_b$residuals)    
  
  expect_equal(fit6_a$coef, fit6_b$coef)
  expect_equal(fit6_a$residuals, fit6_b$residuals)    
  
})

### BIC calculation via 'BICarmaCpp'

test_that("BIC is calculated as intended via BICarmaCpp", {
  
  expect_equal(BIC1_a, BIC1_b)
  expect_equal(BIC2_a, BIC2_b)
  expect_equal(BIC3_a, BIC3_b)
  expect_equal(BIC4_a, BIC4_b)
  expect_equal(BIC5_a, BIC5_b)
  expect_equal(BIC6_a, BIC6_b)  
  
})

### Order selection via 'selectOrderBIC'

BIC_all_b <- matrix(NA, ncol = 4, nrow = 4)
for (p in 0:3) {
  for (q in 0:3) {
    est <- stats::arima(Xt, order = c(p, 0, q))
    BIC_all_b[p + 1, q + 1] <- -2 * est$loglik + ln * (p + q)
  }
}

pq_opt_b <- c(which(BIC_all_b == min(BIC_all_b), arr.ind = TRUE)) - 1
pq_opt_b <- c(pq_opt_b[[1]], 0, pq_opt_b[[2]])
pq_opt_a <- selectOrderBIC(Xt, pmin = 0, pmax = 3, qmin = 0, qmax = 3,
                           armamean = TRUE)

test_that("ARMA model orders are selected correctly following the BIC", {
  
  expect_equal(pq_opt_a, pq_opt_b)
  
})

### Optimal model fitting via 'armaoptCpp'

mopt_b <- stats::arima(Xt, order = pq_opt_b)
mopt_a <- armaoptCpp(Xt, 0, 3, 0, 3, TRUE)

test_that("Selected optimal model return the same fitting results", {
  
  expect_equal(mopt_a$coef, mopt_b$coef)
  expect_equal(mopt_a$residuals, mopt_b$residuals)  
  
})

### Sum of autocovariances is calculated correctly

sum_b <- mopt_b$sigma2 * ((1 + mopt_b$coef[[2]]) / (1 - mopt_b$coef[[1]]))^2
sum_a <- arma_sumacovCpp(mopt_b)

test_that("Sum of ARMA autocovariances is calculated correctly", {
  
  expect_equal(sum_a, sum_b)
  
})