Nothing
tests/testthat/test-varma.R
In ldt: Automated Uncertainty Analysis
test_that("estim.sur estimation works with simulated data", {
set.seed(340)
num_ar <- 1L
num_ma <- 2L
data <- sim.varma(2L, arList = num_ar, maList = num_ma, exoCoef = 2L, nObs = 1000)
res = estim.varma(data = get.data(cbind(data$y,data$x), endogenous = ncol(data$y),
addIntercept = FALSE),
params = c(num_ar, 0, num_ma, 0, 0, 0))
C <- coef(res)
expect_equal(as.numeric(t(coef(res)[c(3,4),])), as.numeric(data$exoCoef), tolerance = 1e-1)
# test changing indexation:
data$x <- data$x[,c(2,1)]
res1 = estim.varma(data = get.data(cbind(data$y,data$x), endogenous = ncol(data$y),
addIntercept = FALSE),
params = c(num_ar, 0, num_ma, 0, 0, 0))
#expect_equal(as.numeric(t(coef(res1)[c(3,4),])), as.numeric(data$exoCoef), tolerance = 1e-1)
# the test does not pass. There are some similar problems in other test.
# My best guess is that it is related to identification when exogenous data presents
})
x <- matrix(c(32.446,44.145,17.062,65.818,76.19,40.408,78.131,
26.695,21.992,68.033,98.872,61.154,71.842,66.922,
31.142,58.429,45.123,80.99,26.345,50.096,36.478,
29.377,27.141,65.037,72.621,63.391,68.125,60.369,
76.384,5.449,99.204,6.87,2.514,98.799,95.082,
6.048,59.287,48.889,21.464,54.972,14.997,27.161,
92.522,65.383,51.852,71.011,55.434,89.082,59.556,
29.524,21.193,2.684,35.457,69.849,76.352,49.455,
18.762,8.492,95.032,39.042,32.517,13.667,91.408,
23.432,56.526,33.531,10.67,72.891,11.796,31.202,
96.893,2.552,82.001,87.786,96.292,93.249,11.688,
19.522,37.55,55.967,97.026,14.017,19.869,60.988,
91.525,33.192,50.666,97.465,58.493,17.033,76.138,
3.432,58.561,69.172,56.453,46.325,63.116,84.577,
12.265,77.277,9.141,69.192,65.464,29.827,8.261,
26.696,94.1,64.958,68.924,97.838,91.389,76.779,
56.448,14.524,33.549,39.059,94.886,98.52,80.476,
2.754,93.605,17.733,37.658,97.567,2.705,74.385,
59.03,10.732,82.043,92.891,69.384,86.848,40.02,
62.295,18.609,61.597,22.438,67.702,83.393,96.283,
64.895,34.39,42.212,52.377,24.745,42.534,64.688,
7.392,82.462,22.022,68.858,55.901,98.156,96.029),
nrow =22, ncol=7)
colnames(x) = paste0("V", c(1:ncol(x)))
test_that("VAR estimation works", {
res = estim.varma(data = get.data(x[,c(1,2)], endogenous = 2), params = c(2,0,0,0,0,0))
#resR = MTS::VAR(x[,1:2],2)
rc1 = 47.01986294876104 # resR$coef[[1]]
rsg = c(498.5024107703099, -56.4870364060901,
-56.4870364060901, 754.8641908858656) # as.numeric(resR$Sigma)
expect_equal(res$estimations$coefs[5,1], rc1, tolerance = 1e-8)
expect_equal(as.numeric(res$estimations$sigma), rsg, tolerance = 1e-8)
#change indexes
res0 = estim.varma(data = get.data(x[,c(2,1)], endogenous = 2), params = c(2,0,0,0,0,0))
expect_equal(res$estimations$coefs[1,1], res0$estimations$coefs[2,2], tolerance = 1e-8) # ??1e-16 fails
expect_equal(res$estimations$sigma[1,1], res0$estimations$sigma[2,2], tolerance = 1e-14)
expect_equal(res$metrics[2,1], res0$metrics[2,1], tolerance = 1e-16)
# test exogenous and endogenous functions:
X <- exogenous(res)
Y <- endogenous(res)
beta <- solve(t(X) %*% X) %*% t(X) %*% Y
expect_equal(res$estimations$coefs, beta)
})
test_that("VAR estimation works with exogenous", {
res = estim.varma(data = get.data(x[,c(1:5)], endogenous = 2), params = c(2,0,0,0,0,0))
X <- exogenous(res)
Y <- endogenous(res)
beta1 <- solve(t(X) %*% X) %*% t(X) %*% Y
expect_equal(res$estimations$coefs, beta1)
#change indexes
res = estim.varma(data = get.data(x[,c(2,1,4,3,5)], endogenous = 2), params = c(2,0,0,0,0,0))
X <- exogenous(res)
Y <- endogenous(res)
beta2 <- solve(t(X) %*% X) %*% t(X) %*% Y
expect_equal(res$estimations$coefs, beta2)
beta1 <- beta1[order(rownames(beta1)),order(colnames(beta1))]
beta2 <- beta2[order(rownames(beta2)),order(colnames(beta2))]
expect_equal(beta1, beta2)
})
test_that("VAR forecast works", {
res = estim.varma(data = get.data(x[,c(1,2)], endogenous = 2), params = c(2,0,0,0,0,0), maxHorizon = 3)
#resR = MTS::VAR(x[,1:2],2)
#forR <- MTS::VARpred(resR, 3);
prd = c(59.71207367716434, 27.85724814539860, 53.92058456057188,
40.35663646224699, 49.31001542183052, 65.52994498334375) #as.numeric(t(forR$pred))
prs = c(22.32716754920583, 27.47479191706219, 22.37005566097166,
28.55614239805317, 22.59875116704332, 32.08105495664848) #as.numeric(t(forR$se.err))
prediction <- predict(res)
expect_equal(as.numeric(t(prediction$means)[,1:3]), prd, tolerance = 1e-8)
#variance
expect_equal(as.numeric(sqrt(t(prediction$vars)[,1:3])), prs, tolerance = 1e-8)
})
test_that("VAR forecast works with NA", {
y = x[,1:2]
y = rbind(c(NA,2),y, c(2,NA))
res = estim.varma(data = get.data(y, endogenous = 2),
params = c(2,0,0,0,0,0), maxHorizon = 3)
#resR = MTS::VAR(x[,1:2],2)
#forR <- MTS::VARpred(resR, 3)
ms1 = c(59.71207367716434, 27.85724814539860, 53.92058456057188,
40.35663646224699, 49.31001542183052, 65.52994498334375) # as.numeric(t(forR$pred))
vs1 = c(22.32716754920583, 27.47479191706219, 22.37005566097166,
28.55614239805317, 22.59875116704332, 32.08105495664848) #as.numeric(t(forR$se.err))
prediction <- predict(res)
expect_equal(as.numeric(t(prediction$means[1:3,])), ms1, tolerance = 1e-8)
#variance
expect_equal(as.numeric(sqrt(t(prediction$vars[1:3,]))), vs1, tolerance = 1e-8)
})
test_that("MA estimation works", {
res = estim.varma(data = get.data(x[,1,drop=FALSE], endogenous = 1),
params = c(0,0,2,0,0,0))
resR = arima(x[,1, drop = FALSE], c(0,0,2), method = "CSS", optim.method = "L-BFGS-B",
transform.pars = FALSE, include.mean = TRUE)
expect_equal(res$estimations$coefs[2:3], as.numeric(resR$coef[1:2]), tolerance = 1e-4)
expect_equal(res$estimations$sigma[[1]], resR$sigma2, tolerance = 1e-8)
})
test_that("MA forecast works", {
res = estim.varma(data = get.data(x[,1,drop=FALSE], endogenous = 1),
params = c(0,0,2,0,0,0), maxHorizon = 3)
resR = arima(x[,1, drop = FALSE], c(0,0,2), method = "CSS", optim.method = "L-BFGS-B",
transform.pars = FALSE, include.mean = TRUE)
forR <- predict(resR, n.ahead = 3)
prediction <- predict(res, actualCount = 0)
expect_equal(as.numeric(forR$pred), as.numeric(prediction$means), tolerance = 1e-6)
#variance
expect_equal(as.numeric(forR$se), sqrt(as.numeric(prediction$vars)), tolerance = 1e-6)
})
test_that("ARMA forecast works", {
res = estim.varma(data = get.data(x[,2,drop=FALSE], endogenous = 1),
params = c(1,0,1,0,0,0), maxHorizon = 3)
resR = arima(x[,2, drop = FALSE], c(1,0,1), method = "CSS", optim.method = "L-BFGS-B",
transform.pars = FALSE, include.mean = TRUE)
forR <- predict(resR, n.ahead = 3)
prediction <- predict(res, actualCount = 0)
expect_equal(as.numeric(forR$pred), as.numeric(prediction$means), tolerance = 1e-3)
#variance
expect_equal(as.numeric(forR$se), sqrt(as.numeric(prediction$vars)), tolerance = 1e-3)
})
test_that("ARIMA forecast works", {
res = estim.varma(data = get.data(x[,5,drop=FALSE], endogenous = 1, addIntercept = FALSE),
params = c(1,1,1,0,0,0), maxHorizon = 3)
resR = arima(x[,5, drop = FALSE], c(1,1,1), method = "CSS", optim.method = "L-BFGS-B",
transform.pars = TRUE, include.mean = TRUE)
forR <- predict(resR, n.ahead = 3)
prediction <- predict(res, actualCount = 0)
expect_equal(as.numeric(forR$pred), as.numeric(prediction$means), tolerance = 1e-3)
#variance
expect_equal(as.numeric(forR$se), sqrt(as.numeric(prediction$vars)), tolerance = 1e-3)
})
test_that("VAR forecast works with PCA for endogenous", {
y = as.matrix(as.data.frame(cbind(as.matrix(x[,1:2]), prcomp(x[,3:ncol(x)], scale. = TRUE)$x))) # orthogonal
pcaOp = get.options.pca()
pcaOp$ignoreFirst = 2
pcaOp$exactCount = 1
res1 = estim.varma(data = get.data(y[,1:3], endogenous = 3),
params = c(1,0,0,0,0,0), maxHorizon = 3)
res2 = estim.varma(data = get.data(x, endogenous = ncol(x)),
params = c(1,0,0,0,0,0), maxHorizon = 3, pcaOptionsY = pcaOp)
expect_equal(res1$gamma, res2$gamma, tolerance = 1e-8)
expect_equal(as.numeric(res1$prediction$means), as.numeric(res2$prediction$means), tolerance = 1e-8)
expect_equal(as.numeric(res1$prediction$vars), as.numeric(res2$prediction$vars), tolerance = 1e-8)
})
test_that("VAR forecast works with PCA for exogenous", {
p=prcomp(x[,3:4], scale. = TRUE)
Z = as.matrix(p$x[,1])
colnames(Z) <- "Z"
newZ = matrix(c(10,11,12, 13,14,15),3,2)
colnames(newZ) = colnames(x[,3:4])
newZp = predict(p,newdata = newZ)
pcaOp = get.options.pca()
pcaOp$ignoreFirst = 0
pcaOp$exactCount = 1
D <- cbind(x[,1:2], Z)
newData = matrix(newZp[,1],3,1)
colnames(newData) <- "Z"
res1 = estim.varma(data = get.data(D, endogenous = 2, newData = newData),
params = c(1,0,0,0,0,0), maxHorizon = 3)
res2 = estim.varma(data = get.data(x[,1:4], endogenous = 2, newData = newZ),
params = c(1,0,0,0,0,0), maxHorizon = 3, pcaOptionsX = pcaOp)
expect_equal(res1$gamma, res2$gamma, tolerance = 1e-8)
expect_equal(res1$prediction$means, res2$prediction$means, tolerance = 1e-8)
expect_equal(res1$prediction$vars, res2$prediction$vars, tolerance = 1e-8)
})
test_that("VAR simulation works", {
Z = x[,1:2]
res = estim.varma(data = get.data(Z, endogenous = 2),
params = c(1,0,0,0,0,0), maxHorizon = 2, simFixSize = 2)
T=nrow(Z)
f1 = estim.varma(data = get.data(Z[1:(T-1),], endogenous = 2),
params = c(1,0,0,0,0,0), maxHorizon = 1)
e1=(abs(f1$prediction$means[,2] - Z[T,])/Z[T,])^2
f2 = estim.varma(data = get.data(Z[1:(T-2),], endogenous = 2),
params = c(1,0,0,0,0,0), maxHorizon = 2)
e2=(abs(f2$prediction$means[,2] - Z[T-1,])/Z[T-1,])^2
e3=(abs(f2$prediction$means[,3] - Z[T,])/Z[T,])^2
expect_equal(as.numeric(sqrt((e1+e2+e3)/3))*100, as.numeric(res$metrics[9,]), tolerance = 1e-14)
# change indexes
Z = x[,c(2,1)]
res0 = estim.varma(data = get.data(Z, endogenous = 2),
params = c(1,0,0,0,0,0), maxHorizon = 2, simFixSize = 2)
expect_equal(res$metrics[,1], res0$metrics[,2], tolerance = 1e-13)
})
test_that("VARMA simulation works", {
Z = x[,1:2]
res = estim.varma(data = get.data(Z, endogenous = 2),
params = c(1,0,1,0,0,0), maxHorizon = 3, simFixSize = 2, simUsePreviousEstim = FALSE)
T=nrow(Z)
f1 = estim.varma(data = get.data(Z[1:(T-1),], endogenous = 2),
params = c(1,0,1,0,0,0), maxHorizon = 1)
e1=(abs(f1$prediction$means[,2] - Z[T,])/Z[T,])^2
f2 = estim.varma(data = get.data(Z[1:(T-2),], endogenous = 2),
params = c(1,0,1,0,0,0), maxHorizon = 2)
e2=(abs(f2$prediction$means[,2] - Z[T-1,])/Z[T-1,])^2
e3=(abs(f2$prediction$means[,3] - Z[T,])/Z[T,])^2
expect_equal(as.numeric(sqrt((e1+e2+e3)/3))*100, as.numeric(res$metrics[9,]), tolerance = 1e-5)
# low tolerance for 'OpenBlas' (TODO: check it)
# change indexes
Z = x[,c(2,1)]
res0 = estim.varma(data = get.data(Z, endogenous = 2),
params = c(1,0,1,0,0,0), maxHorizon = 3, simFixSize = 2, simUsePreviousEstim = FALSE)
expect_equal(res$metrics[,2], res0$metrics[,1], tolerance = 1e-7)
})
test_that("VARMA simulation with lambda in simulation works", {
Z = x[,1:2]
res1 = estim.varma(data = get.data(Z, endogenous = 2, lambdas = c(1, 1)),
params = c(1,0,1,0,0,0), maxHorizon = 3,
simFixSize = 2, simUsePreviousEstim = FALSE)
res2 = estim.varma(data = get.data(Z, endogenous = 2),
params = c(1,0,1,0,0,0), maxHorizon = 3,
simFixSize = 2, simUsePreviousEstim = FALSE)
expect_equal(res1$simulation$validCounts, res2$simulation$validCounts)
# see the note in a similar test in sur test file
})
test_that("ARMA search works for In-Sample", {
skip_on_cran()
res = search.varma(data = get.data(x[,1, drop = FALSE], endogenous = 1, addIntercept = FALSE),
combinations = get.combinations(sizes = c(1),
innerGroups = NULL),
items = get.search.items(bestK = 3, all = TRUE),
maxParams = c(2,2,2,0,0,0),
metrics = get.search.metrics(c("aic")))
sumRes = summary(res, test = TRUE)
expect_equal(sumRes$counts, res$counts)
})
test_that("ARMA search works for In-Sample with exogenous", {
skip_on_cran()
res = search.varma(data = get.data(x[,1:5], endogenous = 1, newData = x),
combinations = get.combinations(sizes = c(1),
innerGroups = list(c(1), c(1,2))),
items = get.search.items(bestK = 3, all = TRUE),
maxParams = c(2,1,2,0,0,0),
metrics = get.search.metrics(c("aic"),c()),
modelChecks = get.search.modelchecks(prediction = FALSE))
sumRes = summary(res, test = TRUE)
expect_equal(sumRes$counts, res$counts)
})
test_that("VARMA search works for In-Sample with exogenous", {
skip_on_cran()
res = search.varma(data = get.data(x[,1:7], endogenous = 3, newData = x),
combinations = get.combinations(sizes = c(1,2),
numTargets = 3,
innerGroups = list(c(1,2))),
maxParams = c(2,2,2,0,0,0),
metrics = get.search.metrics(c("aic", "sic"),c()),
items = get.search.items(all = TRUE),
modelChecks = get.search.modelchecks(prediction = FALSE))
sumRes = summary(res, test = TRUE)
expect_equal(sumRes$counts, res$counts)
})
test_that("VARMA search works when changing Indexes NO exogenous", {
skip_on_cran()
res = search.varma(data = get.data(x[,1:3], endogenous = 3, newData = x),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = NULL),
maxParams = c(2,2,2,0,0,0),
metrics = get.search.metrics(c("aic", "sic"),c()),
items = get.search.items(all = TRUE))
allMetrics = sort(sapply(res$results, function(k){k$value$metric}))
# change place 1 and 2 (both targets)
res0 = search.varma(data = get.data(x[,c(2,1,3)], endogenous = 3),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = NULL),
maxParams = c(2,2,2,0,0,0),
metrics = get.search.metrics(c("aic", "sic"),c()),
items = get.search.items(all = TRUE))
allMetrics0 = sort(sapply(res0$results, function(k){k$value$metric}))
expect_equal(max(abs(allMetrics - allMetrics0)), 0, tolerance = 1e-7)
})
test_that("VARMA search works when changing Indexes WITH exogenous", {
skip_on_cran()
res = search.varma(data = get.data(x[,1:5], endogenous = 3, newData = x),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(1),c(2))),
maxParams = c(0,2,3,0,0,0),
metrics = get.search.metrics(c("sic"),c()),
items = get.search.items(all = TRUE))
allMetrics = sort(sapply(res$results, function(k){k$value$metric}))
res0 = search.varma(data = get.data(x[,c(2,1,3,4,5)], endogenous = 3, newData = x),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(1),c(2))),
maxParams = c(0,2,3,0,0,0),
metrics = get.search.metrics(c("sic"),c()),
items = get.search.items(all = TRUE))
allMetrics0 = sort(sapply(res0$results, function(k){k$value$metric}))
# expect_equal(max(abs(allMetrics - allMetrics0)), 0, tolerance = 1e-7)
# this test fails and need further investigation
# is it related to identification in MA estimation?
})
test_that("V-ARMA search works for Out-Sample", {
skip_on_cran()
res = search.varma(data = get.data(x[,c(1:5)], endogenous = 3, newData = x[5:10,]),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(1,2))),
maxParams = c(2,1,2,0,0,0),
simUsePreviousEstim = FALSE,
maxHorizon = 2,
metrics = get.search.metrics(c(), c("crps", "mae", "rmse"),
horizons = c(1:2), simFixSize = 2),
items = get.search.items(all = TRUE),
modelChecks = get.search.modelchecks(estimation = TRUE,
prediction = TRUE, predictionBound = 5))
sumRes = summary(res, test = TRUE)
expect_equal(sumRes$counts, res$counts)
})
test_that("V-ARMA search works for Out-Sample, no intercept", {
skip_on_cran()
res = search.varma(data = get.data(x[,c(1:7)], endogenous = 3, newData = x[5:10,]),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(3))), #TODO: c(2) fails, why?!
maxParams = c(1,1,1,0,0,0),
simUsePreviousEstim = FALSE,
maxHorizon = 2,
metrics = get.search.metrics(c(), c("mae"),
horizons = c(1:2), simFixSize = 2),
items = get.search.items(all = TRUE),
modelChecks = get.search.modelchecks(estimation = TRUE))
sumRes = summary(res, test = TRUE)
expect_equal(sumRes$counts, res$counts)
})
test_that("VARMA search works when parallel", {
skip_on_cran()
res = search.varma(data = get.data(x[,1:5], endogenous = 3, newData = x),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(1),c(2))),
maxParams = c(0,2,3,0,0,0),
options = get.search.options(parallel = TRUE),
metrics = get.search.metrics(c("sic"),c("mae")),
items = get.search.items(all = TRUE))
allMetrics = sort(sapply(res$results, function(k){k$value$metric}))
res0 = search.varma(data = get.data(x[,c(1:5)], endogenous = 3, newData = x),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(1),c(2))),
maxParams = c(0,2,3,0,0,0),
options = get.search.options(parallel = FALSE),
metrics = get.search.metrics(c("sic"),c("mae")),
items = get.search.items(all = TRUE))
allMetrics0 = sort(sapply(res0$results, function(k){k$value$metric}))
expect_equal(max(abs(allMetrics - allMetrics0)), 0)
})
test_that("VARMA search works with restricted aic", {
skip_on_cran()
res = search.varma(data = get.data(x[,1:7], endogenous = 3, newData = x),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(1,2))),
maxParams = c(0,2,3,0,0,0),
options = get.search.options(parallel = TRUE),
metrics = get.search.metrics(c("sic"),c("mae")),
modelChecks = get.search.modelchecks(maxAic = 220),
items = get.search.items(all = TRUE))
sumRes <- summary(res, test = TRUE)
for (m in sumRes$results){
aic <- as.numeric(m$value$metrics[rownames(m$value$metrics) == "aic",1])
expect_true(aic <= 220)
}
})
test_that("VARMA search works with inclusion weights", {
skip_on_cran()
res = search.varma(data = get.data(x[,c(1:7)], endogenous = 3, newData = x),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(1,2))),
maxParams = c(2,1,2,0,0,0),
metrics = get.search.metrics(c("sic"),c("rmse")),
items = get.search.items(all = TRUE, inclusion = TRUE))
sumRes = summary(res, test = TRUE)
# test fails in some cases when 3 or 4 are in innerGroups = list(c(1,2)) ?!!
#allMetrics = sapply(res$results, function(k){k$value$metric})
#which(abs(allMetrics -3.92371108382458)<1e-6 )
inclusion = matrix(0,8,2, dimnames = list(colnames(res$info$data$data), NULL))
for (m in res$results[which(sapply(res$results,
function(r) r$evalName == "sic" && r$typeName == "model" && r$targetName == "V1"))]){
for (d in (m$value$endogenous)){
inclusion[d,1] = inclusion[d,1] + m$value$weight
inclusion[d,2] = inclusion[d,2] + 1
}
for (d in (m$value$exogenous)){
inclusion[d,1] = inclusion[d,1] + m$value$weight
inclusion[d,2] = inclusion[d,2] + 1
}
}
inclusion[,1] = inclusion[,1]/inclusion[,2]
searchInclusion = res$results[which(sapply(res$results,
function(r) r$evalName == "sic" && r$targetName == "V1" && r$typeName == "inclusion"))]
expect_equal(as.numeric(searchInclusion[[1]]$value), as.numeric(inclusion), tolerance = 1e-10)
})
test_that("VARMA search works with predictions (bests)", {
skip_on_cran()
res = search.varma(data = get.data(x[,1:7], endogenous = 3, newData = x[c(8,9,10),]),
combinations = get.combinations(sizes = c(1,2,3),
numTargets = 2,
innerGroups = NULL),
maxParams = c(2,2,2,0,0,0),
maxHorizon = 3,
simUsePreviousEstim = FALSE,
options = get.search.options(parallel = FALSE),
metrics = get.search.metrics(c("sic"),c("mae"),horizons = c(1:3)),
modelChecks = get.search.modelchecks(predictionBound = 300),
items = get.search.items(all = TRUE, type1 = TRUE, bestK = 3))
sumRes = summary(res, test = TRUE)
expect_equal(sumRes$counts, res$counts)
})
test_that("VARMA search works with predictions (cdfs)", {
skip_on_cran()
res = search.varma(data = get.data(x[,1:7], endogenous = 3, newData = x[c(8,9,10),]),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(1,2))),
maxParams = c(2,1,2,0,0,0),
maxHorizon = 3,
simUsePreviousEstim = FALSE,
metrics = get.search.metrics(c("sic"),c("rmse"),
horizons = c(1:3), simFixSize = 2),
items = get.search.items(all = TRUE, type1 = TRUE, cdfs = c(0,1,0)),
modelChecks = get.search.modelchecks(estimation = FALSE, prediction = FALSE, predictionBound = 0))
sumRes <- summary(res, test = TRUE)
h = 2
sum = 0
c = 0
cc=0
i = 0
for (m in sumRes$results){
i = i + 1
if (m$evalName != "rmse" || m$typeName != "model" || m$targetName != "V1")
next()
hh= h + m$value$prediction$startIndex - 1
coef = m$value$prediction$means["V1",hh]
#print(coef)
sd = sqrt(m$value$prediction$vars["V1",hh])
w = res$results[[i]]$value$weight
sum = sum + w * pnorm(0,coef,sd) # note the NORMAL dist.
c=c+w
cc=cc+1
}
cdfs = res$results[which(sapply(res$results,
function(r) r$evalName == "rmse" && r$targetName == "V1" && r$typeName == "cdf"))]
expect_equal(cdfs[[1]]$value[2,1], sum/c, tolerance = 1e-10)
expect_equal(cdfs[[1]]$value[2,3], c, tolerance = 1e-10)
# test does not pass the second column for count (actual:48, res[2,2]: 46) in Debian
# note that CDF calculations are weight based
# same in mixture4
})
test_that("VARMA search works with predictions (extreme bounds)", {
skip_on_cran()
res = search.varma(data = get.data(x[,1:7], endogenous = 3, newData = x[c(8,9,10),]),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(1,2))),
maxParams = c(2,1,2,0,0,0),
maxHorizon = 3,
simUsePreviousEstim = FALSE,
metrics = get.search.metrics(c("sic"),c("rmse"),
horizons = c(1:3), simFixSize = 2),
items = get.search.items(all = TRUE, type1 = TRUE, extremeMultiplier = 2),
modelChecks = get.search.modelchecks(estimation = FALSE, prediction = FALSE, predictionBound = 0))
sumRes <- summary(res, test = TRUE)
h = 2
mn = Inf
mx = -Inf
for (m in sumRes$results){
if (m$evalName != "rmse" || m$typeName != "model" || m$targetName != "V1")
next()
hh= h + m$value$prediction$startIndex - 1
coef = m$value$prediction$means["V1",hh]
sd = sqrt(m$value$prediction$vars["V1",hh])
mn = min(mn,coef-2*sd)
mx = max(mx,coef+2*sd)
}
extremeB = res$results[which(sapply(res$results,
function(r) r$evalName == "rmse" && r$targetName == "V1" && r$typeName == "extreme bound"))]
expect_equal(extremeB[[1]]$value[2,1], mn, tolerance = 1e-10)
expect_equal(extremeB[[1]]$value[2,2], mx, tolerance = 1e-10)
})
test_that("VARMA search works with predictions (mixture)", {
skip_on_cran()
res = search.varma(data = get.data(x[,1:7], endogenous = 3, newData = x[c(8,9,10),]),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(1,2))),
maxParams = c(2,1,2,0,0,0),
maxHorizon = 3,
simUsePreviousEstim = FALSE,
metrics = get.search.metrics(c("sic"),c("rmse"),
horizons = c(1:3), simFixSize = 2),
items = get.search.items(all = TRUE, type1 = TRUE, mixture4 = TRUE),
modelChecks = get.search.modelchecks(estimation = FALSE, prediction = FALSE, predictionBound = 0))
sumRes <- summary(res, test = TRUE)
h = 2
coefs = c()
vars = c()
weights = c()
i <- 0
for (m in sumRes$results){
i <- i + 1
if (m$evalName != "rmse" || m$typeName != "model" || m$targetName != "V1")
next()
hh= h + m$value$prediction$startIndex - 1
coefs = append(coefs,m$value$prediction$means["V1",hh])
vars = append(vars, m$value$prediction$vars["V1",hh])
weights = append(weights, res$results[[i]]$value$weight)
}
mixture = res$results[which(sapply(res$results,
function(r) r$evalName == "rmse" && r$targetName == "V1" && r$typeName == "mixture"))]
# note that we need weighted mean, variance, etc. assuming normal distribution
me = weighted.mean(coefs, weights)
expect_equal(mixture[[1]]$value[2,1], me, tolerance = 1e-14)
# TODO : compare weighted variance, skewness, kurtosis assuming normality
# of course, its better to .Call the running statistics, test it, and use it here
#len = length(coefs)
# expect_equal(mixture[[1]]$value[2,5], len)
# test does not pass the second column for count (actual:48, res[2,2]: 46) in Debian
# note that mixture calculations are weight based
# same in CDF
})
test_that("estim.varma SplitSearch works (no subsetting)", {
skip_on_cran()
data = data = get.data(x[,1:7], endogenous = 3, newData = x[c(8,9,10),])
combinations = get.combinations(numTargets = 3, innerGroups = list(c(1), c(1,2), c(1,3)))
items = get.search.items(inclusion = TRUE
#, all = TRUE
, bestK = 4
, type1 = TRUE
, cdfs = c(0,0.3)
, mixture4 = TRUE
, extremeMultiplier = 2
)
metrics = get.search.metrics(c("sic", "aic"),
horizons = c(1:3), simFixSize = 2) # don't test with out-of-sample metrics. It seems we have different model with equal weights (the result change by repeating the call ?!)
options = get.search.options(FALSE,
#reportInterval = 1
)
combinations$sizes <- c(1, 2, 3)
whole = search.varma(data = data,
combinations = combinations,
maxParams = c(2,1,2,0,0,0),
maxHorizon = 3,
simUsePreviousEstim = FALSE,
items = items,
metrics = metrics,
options = options)
combinations$sizes <- list(c(1, 2), c(3))
combinations$stepsNumVariables <- c(NA, NA)
split = search.varma(data = data,
combinations = combinations,
maxParams = c(2,1,2,0,0,0),
maxHorizon = 3,
simUsePreviousEstim = FALSE,
items = items,
metrics = metrics,
options = options)
expect_equal(whole$counts, split$counts)
expect_equal(length(whole$results), length(split$results))
pastedList_w <- unlist(lapply(whole$results, function(x) paste(x[1:4], collapse = "")))
pastedList_s <- unlist(lapply(split$results, function(x) paste(x[1:4], collapse = "")))
sortedList_w <- whole$results[order(pastedList_w)]
sortedList_s <- split$results[order(pastedList_s)]
for (i in 1:length(sortedList_w)){
if (sortedList_s[[i]]$typeName == "mixture"){
expect_equal(sortedList_s[[i]]$value[,c(1:3,5,6)], sortedList_w[[i]]$value[,c(1:3,5,6)])
expect_equal(sortedList_s[[i]]$value[,c(4)], sortedList_w[[i]]$value[,c(4)], tolerance = 0.1)
}
else
expect_equal(sortedList_s[[i]]$value, sortedList_w[[i]]$value)
}
})
Try the ldt package in your browser
Any scripts or data that you put into this service are public.
ldt documentation built on Sept. 11, 2024, 5:25 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
test_that("estim.sur estimation works with simulated data", {
set.seed(340)
num_ar <- 1L
num_ma <- 2L
data <- sim.varma(2L, arList = num_ar, maList = num_ma, exoCoef = 2L, nObs = 1000)
res = estim.varma(data = get.data(cbind(data$y,data$x), endogenous = ncol(data$y),
addIntercept = FALSE),
params = c(num_ar, 0, num_ma, 0, 0, 0))
C <- coef(res)
expect_equal(as.numeric(t(coef(res)[c(3,4),])), as.numeric(data$exoCoef), tolerance = 1e-1)
# test changing indexation:
data$x <- data$x[,c(2,1)]
res1 = estim.varma(data = get.data(cbind(data$y,data$x), endogenous = ncol(data$y),
addIntercept = FALSE),
params = c(num_ar, 0, num_ma, 0, 0, 0))
#expect_equal(as.numeric(t(coef(res1)[c(3,4),])), as.numeric(data$exoCoef), tolerance = 1e-1)
# the test does not pass. There are some similar problems in other test.
# My best guess is that it is related to identification when exogenous data presents
})
x <- matrix(c(32.446,44.145,17.062,65.818,76.19,40.408,78.131,
26.695,21.992,68.033,98.872,61.154,71.842,66.922,
31.142,58.429,45.123,80.99,26.345,50.096,36.478,
29.377,27.141,65.037,72.621,63.391,68.125,60.369,
76.384,5.449,99.204,6.87,2.514,98.799,95.082,
6.048,59.287,48.889,21.464,54.972,14.997,27.161,
92.522,65.383,51.852,71.011,55.434,89.082,59.556,
29.524,21.193,2.684,35.457,69.849,76.352,49.455,
18.762,8.492,95.032,39.042,32.517,13.667,91.408,
23.432,56.526,33.531,10.67,72.891,11.796,31.202,
96.893,2.552,82.001,87.786,96.292,93.249,11.688,
19.522,37.55,55.967,97.026,14.017,19.869,60.988,
91.525,33.192,50.666,97.465,58.493,17.033,76.138,
3.432,58.561,69.172,56.453,46.325,63.116,84.577,
12.265,77.277,9.141,69.192,65.464,29.827,8.261,
26.696,94.1,64.958,68.924,97.838,91.389,76.779,
56.448,14.524,33.549,39.059,94.886,98.52,80.476,
2.754,93.605,17.733,37.658,97.567,2.705,74.385,
59.03,10.732,82.043,92.891,69.384,86.848,40.02,
62.295,18.609,61.597,22.438,67.702,83.393,96.283,
64.895,34.39,42.212,52.377,24.745,42.534,64.688,
7.392,82.462,22.022,68.858,55.901,98.156,96.029),
nrow =22, ncol=7)
colnames(x) = paste0("V", c(1:ncol(x)))
test_that("VAR estimation works", {
res = estim.varma(data = get.data(x[,c(1,2)], endogenous = 2), params = c(2,0,0,0,0,0))
#resR = MTS::VAR(x[,1:2],2)
rc1 = 47.01986294876104 # resR$coef[[1]]
rsg = c(498.5024107703099, -56.4870364060901,
-56.4870364060901, 754.8641908858656) # as.numeric(resR$Sigma)
expect_equal(res$estimations$coefs[5,1], rc1, tolerance = 1e-8)
expect_equal(as.numeric(res$estimations$sigma), rsg, tolerance = 1e-8)
#change indexes
res0 = estim.varma(data = get.data(x[,c(2,1)], endogenous = 2), params = c(2,0,0,0,0,0))
expect_equal(res$estimations$coefs[1,1], res0$estimations$coefs[2,2], tolerance = 1e-8) # ??1e-16 fails
expect_equal(res$estimations$sigma[1,1], res0$estimations$sigma[2,2], tolerance = 1e-14)
expect_equal(res$metrics[2,1], res0$metrics[2,1], tolerance = 1e-16)
# test exogenous and endogenous functions:
X <- exogenous(res)
Y <- endogenous(res)
beta <- solve(t(X) %*% X) %*% t(X) %*% Y
expect_equal(res$estimations$coefs, beta)
})
test_that("VAR estimation works with exogenous", {
res = estim.varma(data = get.data(x[,c(1:5)], endogenous = 2), params = c(2,0,0,0,0,0))
X <- exogenous(res)
Y <- endogenous(res)
beta1 <- solve(t(X) %*% X) %*% t(X) %*% Y
expect_equal(res$estimations$coefs, beta1)
#change indexes
res = estim.varma(data = get.data(x[,c(2,1,4,3,5)], endogenous = 2), params = c(2,0,0,0,0,0))
X <- exogenous(res)
Y <- endogenous(res)
beta2 <- solve(t(X) %*% X) %*% t(X) %*% Y
expect_equal(res$estimations$coefs, beta2)
beta1 <- beta1[order(rownames(beta1)),order(colnames(beta1))]
beta2 <- beta2[order(rownames(beta2)),order(colnames(beta2))]
expect_equal(beta1, beta2)
})
test_that("VAR forecast works", {
res = estim.varma(data = get.data(x[,c(1,2)], endogenous = 2), params = c(2,0,0,0,0,0), maxHorizon = 3)
#resR = MTS::VAR(x[,1:2],2)
#forR <- MTS::VARpred(resR, 3);
prd = c(59.71207367716434, 27.85724814539860, 53.92058456057188,
40.35663646224699, 49.31001542183052, 65.52994498334375) #as.numeric(t(forR$pred))
prs = c(22.32716754920583, 27.47479191706219, 22.37005566097166,
28.55614239805317, 22.59875116704332, 32.08105495664848) #as.numeric(t(forR$se.err))
prediction <- predict(res)
expect_equal(as.numeric(t(prediction$means)[,1:3]), prd, tolerance = 1e-8)
#variance
expect_equal(as.numeric(sqrt(t(prediction$vars)[,1:3])), prs, tolerance = 1e-8)
})
test_that("VAR forecast works with NA", {
y = x[,1:2]
y = rbind(c(NA,2),y, c(2,NA))
res = estim.varma(data = get.data(y, endogenous = 2),
params = c(2,0,0,0,0,0), maxHorizon = 3)
#resR = MTS::VAR(x[,1:2],2)
#forR <- MTS::VARpred(resR, 3)
ms1 = c(59.71207367716434, 27.85724814539860, 53.92058456057188,
40.35663646224699, 49.31001542183052, 65.52994498334375) # as.numeric(t(forR$pred))
vs1 = c(22.32716754920583, 27.47479191706219, 22.37005566097166,
28.55614239805317, 22.59875116704332, 32.08105495664848) #as.numeric(t(forR$se.err))
prediction <- predict(res)
expect_equal(as.numeric(t(prediction$means[1:3,])), ms1, tolerance = 1e-8)
#variance
expect_equal(as.numeric(sqrt(t(prediction$vars[1:3,]))), vs1, tolerance = 1e-8)
})
test_that("MA estimation works", {
res = estim.varma(data = get.data(x[,1,drop=FALSE], endogenous = 1),
params = c(0,0,2,0,0,0))
resR = arima(x[,1, drop = FALSE], c(0,0,2), method = "CSS", optim.method = "L-BFGS-B",
transform.pars = FALSE, include.mean = TRUE)
expect_equal(res$estimations$coefs[2:3], as.numeric(resR$coef[1:2]), tolerance = 1e-4)
expect_equal(res$estimations$sigma[[1]], resR$sigma2, tolerance = 1e-8)
})
test_that("MA forecast works", {
res = estim.varma(data = get.data(x[,1,drop=FALSE], endogenous = 1),
params = c(0,0,2,0,0,0), maxHorizon = 3)
resR = arima(x[,1, drop = FALSE], c(0,0,2), method = "CSS", optim.method = "L-BFGS-B",
transform.pars = FALSE, include.mean = TRUE)
forR <- predict(resR, n.ahead = 3)
prediction <- predict(res, actualCount = 0)
expect_equal(as.numeric(forR$pred), as.numeric(prediction$means), tolerance = 1e-6)
#variance
expect_equal(as.numeric(forR$se), sqrt(as.numeric(prediction$vars)), tolerance = 1e-6)
})
test_that("ARMA forecast works", {
res = estim.varma(data = get.data(x[,2,drop=FALSE], endogenous = 1),
params = c(1,0,1,0,0,0), maxHorizon = 3)
resR = arima(x[,2, drop = FALSE], c(1,0,1), method = "CSS", optim.method = "L-BFGS-B",
transform.pars = FALSE, include.mean = TRUE)
forR <- predict(resR, n.ahead = 3)
prediction <- predict(res, actualCount = 0)
expect_equal(as.numeric(forR$pred), as.numeric(prediction$means), tolerance = 1e-3)
#variance
expect_equal(as.numeric(forR$se), sqrt(as.numeric(prediction$vars)), tolerance = 1e-3)
})
test_that("ARIMA forecast works", {
res = estim.varma(data = get.data(x[,5,drop=FALSE], endogenous = 1, addIntercept = FALSE),
params = c(1,1,1,0,0,0), maxHorizon = 3)
resR = arima(x[,5, drop = FALSE], c(1,1,1), method = "CSS", optim.method = "L-BFGS-B",
transform.pars = TRUE, include.mean = TRUE)
forR <- predict(resR, n.ahead = 3)
prediction <- predict(res, actualCount = 0)
expect_equal(as.numeric(forR$pred), as.numeric(prediction$means), tolerance = 1e-3)
#variance
expect_equal(as.numeric(forR$se), sqrt(as.numeric(prediction$vars)), tolerance = 1e-3)
})
test_that("VAR forecast works with PCA for endogenous", {
y = as.matrix(as.data.frame(cbind(as.matrix(x[,1:2]), prcomp(x[,3:ncol(x)], scale. = TRUE)$x))) # orthogonal
pcaOp = get.options.pca()
pcaOp$ignoreFirst = 2
pcaOp$exactCount = 1
res1 = estim.varma(data = get.data(y[,1:3], endogenous = 3),
params = c(1,0,0,0,0,0), maxHorizon = 3)
res2 = estim.varma(data = get.data(x, endogenous = ncol(x)),
params = c(1,0,0,0,0,0), maxHorizon = 3, pcaOptionsY = pcaOp)
expect_equal(res1$gamma, res2$gamma, tolerance = 1e-8)
expect_equal(as.numeric(res1$prediction$means), as.numeric(res2$prediction$means), tolerance = 1e-8)
expect_equal(as.numeric(res1$prediction$vars), as.numeric(res2$prediction$vars), tolerance = 1e-8)
})
test_that("VAR forecast works with PCA for exogenous", {
p=prcomp(x[,3:4], scale. = TRUE)
Z = as.matrix(p$x[,1])
colnames(Z) <- "Z"
newZ = matrix(c(10,11,12, 13,14,15),3,2)
colnames(newZ) = colnames(x[,3:4])
newZp = predict(p,newdata = newZ)
pcaOp = get.options.pca()
pcaOp$ignoreFirst = 0
pcaOp$exactCount = 1
D <- cbind(x[,1:2], Z)
newData = matrix(newZp[,1],3,1)
colnames(newData) <- "Z"
res1 = estim.varma(data = get.data(D, endogenous = 2, newData = newData),
params = c(1,0,0,0,0,0), maxHorizon = 3)
res2 = estim.varma(data = get.data(x[,1:4], endogenous = 2, newData = newZ),
params = c(1,0,0,0,0,0), maxHorizon = 3, pcaOptionsX = pcaOp)
expect_equal(res1$gamma, res2$gamma, tolerance = 1e-8)
expect_equal(res1$prediction$means, res2$prediction$means, tolerance = 1e-8)
expect_equal(res1$prediction$vars, res2$prediction$vars, tolerance = 1e-8)
})
test_that("VAR simulation works", {
Z = x[,1:2]
res = estim.varma(data = get.data(Z, endogenous = 2),
params = c(1,0,0,0,0,0), maxHorizon = 2, simFixSize = 2)
T=nrow(Z)
f1 = estim.varma(data = get.data(Z[1:(T-1),], endogenous = 2),
params = c(1,0,0,0,0,0), maxHorizon = 1)
e1=(abs(f1$prediction$means[,2] - Z[T,])/Z[T,])^2
f2 = estim.varma(data = get.data(Z[1:(T-2),], endogenous = 2),
params = c(1,0,0,0,0,0), maxHorizon = 2)
e2=(abs(f2$prediction$means[,2] - Z[T-1,])/Z[T-1,])^2
e3=(abs(f2$prediction$means[,3] - Z[T,])/Z[T,])^2
expect_equal(as.numeric(sqrt((e1+e2+e3)/3))*100, as.numeric(res$metrics[9,]), tolerance = 1e-14)
# change indexes
Z = x[,c(2,1)]
res0 = estim.varma(data = get.data(Z, endogenous = 2),
params = c(1,0,0,0,0,0), maxHorizon = 2, simFixSize = 2)
expect_equal(res$metrics[,1], res0$metrics[,2], tolerance = 1e-13)
})
test_that("VARMA simulation works", {
Z = x[,1:2]
res = estim.varma(data = get.data(Z, endogenous = 2),
params = c(1,0,1,0,0,0), maxHorizon = 3, simFixSize = 2, simUsePreviousEstim = FALSE)
T=nrow(Z)
f1 = estim.varma(data = get.data(Z[1:(T-1),], endogenous = 2),
params = c(1,0,1,0,0,0), maxHorizon = 1)
e1=(abs(f1$prediction$means[,2] - Z[T,])/Z[T,])^2
f2 = estim.varma(data = get.data(Z[1:(T-2),], endogenous = 2),
params = c(1,0,1,0,0,0), maxHorizon = 2)
e2=(abs(f2$prediction$means[,2] - Z[T-1,])/Z[T-1,])^2
e3=(abs(f2$prediction$means[,3] - Z[T,])/Z[T,])^2
expect_equal(as.numeric(sqrt((e1+e2+e3)/3))*100, as.numeric(res$metrics[9,]), tolerance = 1e-5)
# low tolerance for 'OpenBlas' (TODO: check it)
# change indexes
Z = x[,c(2,1)]
res0 = estim.varma(data = get.data(Z, endogenous = 2),
params = c(1,0,1,0,0,0), maxHorizon = 3, simFixSize = 2, simUsePreviousEstim = FALSE)
expect_equal(res$metrics[,2], res0$metrics[,1], tolerance = 1e-7)
})
test_that("VARMA simulation with lambda in simulation works", {
Z = x[,1:2]
res1 = estim.varma(data = get.data(Z, endogenous = 2, lambdas = c(1, 1)),
params = c(1,0,1,0,0,0), maxHorizon = 3,
simFixSize = 2, simUsePreviousEstim = FALSE)
res2 = estim.varma(data = get.data(Z, endogenous = 2),
params = c(1,0,1,0,0,0), maxHorizon = 3,
simFixSize = 2, simUsePreviousEstim = FALSE)
expect_equal(res1$simulation$validCounts, res2$simulation$validCounts)
# see the note in a similar test in sur test file
})
test_that("ARMA search works for In-Sample", {
skip_on_cran()
res = search.varma(data = get.data(x[,1, drop = FALSE], endogenous = 1, addIntercept = FALSE),
combinations = get.combinations(sizes = c(1),
innerGroups = NULL),
items = get.search.items(bestK = 3, all = TRUE),
maxParams = c(2,2,2,0,0,0),
metrics = get.search.metrics(c("aic")))
sumRes = summary(res, test = TRUE)
expect_equal(sumRes$counts, res$counts)
})
test_that("ARMA search works for In-Sample with exogenous", {
skip_on_cran()
res = search.varma(data = get.data(x[,1:5], endogenous = 1, newData = x),
combinations = get.combinations(sizes = c(1),
innerGroups = list(c(1), c(1,2))),
items = get.search.items(bestK = 3, all = TRUE),
maxParams = c(2,1,2,0,0,0),
metrics = get.search.metrics(c("aic"),c()),
modelChecks = get.search.modelchecks(prediction = FALSE))
sumRes = summary(res, test = TRUE)
expect_equal(sumRes$counts, res$counts)
})
test_that("VARMA search works for In-Sample with exogenous", {
skip_on_cran()
res = search.varma(data = get.data(x[,1:7], endogenous = 3, newData = x),
combinations = get.combinations(sizes = c(1,2),
numTargets = 3,
innerGroups = list(c(1,2))),
maxParams = c(2,2,2,0,0,0),
metrics = get.search.metrics(c("aic", "sic"),c()),
items = get.search.items(all = TRUE),
modelChecks = get.search.modelchecks(prediction = FALSE))
sumRes = summary(res, test = TRUE)
expect_equal(sumRes$counts, res$counts)
})
test_that("VARMA search works when changing Indexes NO exogenous", {
skip_on_cran()
res = search.varma(data = get.data(x[,1:3], endogenous = 3, newData = x),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = NULL),
maxParams = c(2,2,2,0,0,0),
metrics = get.search.metrics(c("aic", "sic"),c()),
items = get.search.items(all = TRUE))
allMetrics = sort(sapply(res$results, function(k){k$value$metric}))
# change place 1 and 2 (both targets)
res0 = search.varma(data = get.data(x[,c(2,1,3)], endogenous = 3),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = NULL),
maxParams = c(2,2,2,0,0,0),
metrics = get.search.metrics(c("aic", "sic"),c()),
items = get.search.items(all = TRUE))
allMetrics0 = sort(sapply(res0$results, function(k){k$value$metric}))
expect_equal(max(abs(allMetrics - allMetrics0)), 0, tolerance = 1e-7)
})
test_that("VARMA search works when changing Indexes WITH exogenous", {
skip_on_cran()
res = search.varma(data = get.data(x[,1:5], endogenous = 3, newData = x),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(1),c(2))),
maxParams = c(0,2,3,0,0,0),
metrics = get.search.metrics(c("sic"),c()),
items = get.search.items(all = TRUE))
allMetrics = sort(sapply(res$results, function(k){k$value$metric}))
res0 = search.varma(data = get.data(x[,c(2,1,3,4,5)], endogenous = 3, newData = x),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(1),c(2))),
maxParams = c(0,2,3,0,0,0),
metrics = get.search.metrics(c("sic"),c()),
items = get.search.items(all = TRUE))
allMetrics0 = sort(sapply(res0$results, function(k){k$value$metric}))
# expect_equal(max(abs(allMetrics - allMetrics0)), 0, tolerance = 1e-7)
# this test fails and need further investigation
# is it related to identification in MA estimation?
})
test_that("V-ARMA search works for Out-Sample", {
skip_on_cran()
res = search.varma(data = get.data(x[,c(1:5)], endogenous = 3, newData = x[5:10,]),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(1,2))),
maxParams = c(2,1,2,0,0,0),
simUsePreviousEstim = FALSE,
maxHorizon = 2,
metrics = get.search.metrics(c(), c("crps", "mae", "rmse"),
horizons = c(1:2), simFixSize = 2),
items = get.search.items(all = TRUE),
modelChecks = get.search.modelchecks(estimation = TRUE,
prediction = TRUE, predictionBound = 5))
sumRes = summary(res, test = TRUE)
expect_equal(sumRes$counts, res$counts)
})
test_that("V-ARMA search works for Out-Sample, no intercept", {
skip_on_cran()
res = search.varma(data = get.data(x[,c(1:7)], endogenous = 3, newData = x[5:10,]),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(3))), #TODO: c(2) fails, why?!
maxParams = c(1,1,1,0,0,0),
simUsePreviousEstim = FALSE,
maxHorizon = 2,
metrics = get.search.metrics(c(), c("mae"),
horizons = c(1:2), simFixSize = 2),
items = get.search.items(all = TRUE),
modelChecks = get.search.modelchecks(estimation = TRUE))
sumRes = summary(res, test = TRUE)
expect_equal(sumRes$counts, res$counts)
})
test_that("VARMA search works when parallel", {
skip_on_cran()
res = search.varma(data = get.data(x[,1:5], endogenous = 3, newData = x),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(1),c(2))),
maxParams = c(0,2,3,0,0,0),
options = get.search.options(parallel = TRUE),
metrics = get.search.metrics(c("sic"),c("mae")),
items = get.search.items(all = TRUE))
allMetrics = sort(sapply(res$results, function(k){k$value$metric}))
res0 = search.varma(data = get.data(x[,c(1:5)], endogenous = 3, newData = x),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(1),c(2))),
maxParams = c(0,2,3,0,0,0),
options = get.search.options(parallel = FALSE),
metrics = get.search.metrics(c("sic"),c("mae")),
items = get.search.items(all = TRUE))
allMetrics0 = sort(sapply(res0$results, function(k){k$value$metric}))
expect_equal(max(abs(allMetrics - allMetrics0)), 0)
})
test_that("VARMA search works with restricted aic", {
skip_on_cran()
res = search.varma(data = get.data(x[,1:7], endogenous = 3, newData = x),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(1,2))),
maxParams = c(0,2,3,0,0,0),
options = get.search.options(parallel = TRUE),
metrics = get.search.metrics(c("sic"),c("mae")),
modelChecks = get.search.modelchecks(maxAic = 220),
items = get.search.items(all = TRUE))
sumRes <- summary(res, test = TRUE)
for (m in sumRes$results){
aic <- as.numeric(m$value$metrics[rownames(m$value$metrics) == "aic",1])
expect_true(aic <= 220)
}
})
test_that("VARMA search works with inclusion weights", {
skip_on_cran()
res = search.varma(data = get.data(x[,c(1:7)], endogenous = 3, newData = x),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(1,2))),
maxParams = c(2,1,2,0,0,0),
metrics = get.search.metrics(c("sic"),c("rmse")),
items = get.search.items(all = TRUE, inclusion = TRUE))
sumRes = summary(res, test = TRUE)
# test fails in some cases when 3 or 4 are in innerGroups = list(c(1,2)) ?!!
#allMetrics = sapply(res$results, function(k){k$value$metric})
#which(abs(allMetrics -3.92371108382458)<1e-6 )
inclusion = matrix(0,8,2, dimnames = list(colnames(res$info$data$data), NULL))
for (m in res$results[which(sapply(res$results,
function(r) r$evalName == "sic" && r$typeName == "model" && r$targetName == "V1"))]){
for (d in (m$value$endogenous)){
inclusion[d,1] = inclusion[d,1] + m$value$weight
inclusion[d,2] = inclusion[d,2] + 1
}
for (d in (m$value$exogenous)){
inclusion[d,1] = inclusion[d,1] + m$value$weight
inclusion[d,2] = inclusion[d,2] + 1
}
}
inclusion[,1] = inclusion[,1]/inclusion[,2]
searchInclusion = res$results[which(sapply(res$results,
function(r) r$evalName == "sic" && r$targetName == "V1" && r$typeName == "inclusion"))]
expect_equal(as.numeric(searchInclusion[[1]]$value), as.numeric(inclusion), tolerance = 1e-10)
})
test_that("VARMA search works with predictions (bests)", {
skip_on_cran()
res = search.varma(data = get.data(x[,1:7], endogenous = 3, newData = x[c(8,9,10),]),
combinations = get.combinations(sizes = c(1,2,3),
numTargets = 2,
innerGroups = NULL),
maxParams = c(2,2,2,0,0,0),
maxHorizon = 3,
simUsePreviousEstim = FALSE,
options = get.search.options(parallel = FALSE),
metrics = get.search.metrics(c("sic"),c("mae"),horizons = c(1:3)),
modelChecks = get.search.modelchecks(predictionBound = 300),
items = get.search.items(all = TRUE, type1 = TRUE, bestK = 3))
sumRes = summary(res, test = TRUE)
expect_equal(sumRes$counts, res$counts)
})
test_that("VARMA search works with predictions (cdfs)", {
skip_on_cran()
res = search.varma(data = get.data(x[,1:7], endogenous = 3, newData = x[c(8,9,10),]),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(1,2))),
maxParams = c(2,1,2,0,0,0),
maxHorizon = 3,
simUsePreviousEstim = FALSE,
metrics = get.search.metrics(c("sic"),c("rmse"),
horizons = c(1:3), simFixSize = 2),
items = get.search.items(all = TRUE, type1 = TRUE, cdfs = c(0,1,0)),
modelChecks = get.search.modelchecks(estimation = FALSE, prediction = FALSE, predictionBound = 0))
sumRes <- summary(res, test = TRUE)
h = 2
sum = 0
c = 0
cc=0
i = 0
for (m in sumRes$results){
i = i + 1
if (m$evalName != "rmse" || m$typeName != "model" || m$targetName != "V1")
next()
hh= h + m$value$prediction$startIndex - 1
coef = m$value$prediction$means["V1",hh]
#print(coef)
sd = sqrt(m$value$prediction$vars["V1",hh])
w = res$results[[i]]$value$weight
sum = sum + w * pnorm(0,coef,sd) # note the NORMAL dist.
c=c+w
cc=cc+1
}
cdfs = res$results[which(sapply(res$results,
function(r) r$evalName == "rmse" && r$targetName == "V1" && r$typeName == "cdf"))]
expect_equal(cdfs[[1]]$value[2,1], sum/c, tolerance = 1e-10)
expect_equal(cdfs[[1]]$value[2,3], c, tolerance = 1e-10)
# test does not pass the second column for count (actual:48, res[2,2]: 46) in Debian
# note that CDF calculations are weight based
# same in mixture4
})
test_that("VARMA search works with predictions (extreme bounds)", {
skip_on_cran()
res = search.varma(data = get.data(x[,1:7], endogenous = 3, newData = x[c(8,9,10),]),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(1,2))),
maxParams = c(2,1,2,0,0,0),
maxHorizon = 3,
simUsePreviousEstim = FALSE,
metrics = get.search.metrics(c("sic"),c("rmse"),
horizons = c(1:3), simFixSize = 2),
items = get.search.items(all = TRUE, type1 = TRUE, extremeMultiplier = 2),
modelChecks = get.search.modelchecks(estimation = FALSE, prediction = FALSE, predictionBound = 0))
sumRes <- summary(res, test = TRUE)
h = 2
mn = Inf
mx = -Inf
for (m in sumRes$results){
if (m$evalName != "rmse" || m$typeName != "model" || m$targetName != "V1")
next()
hh= h + m$value$prediction$startIndex - 1
coef = m$value$prediction$means["V1",hh]
sd = sqrt(m$value$prediction$vars["V1",hh])
mn = min(mn,coef-2*sd)
mx = max(mx,coef+2*sd)
}
extremeB = res$results[which(sapply(res$results,
function(r) r$evalName == "rmse" && r$targetName == "V1" && r$typeName == "extreme bound"))]
expect_equal(extremeB[[1]]$value[2,1], mn, tolerance = 1e-10)
expect_equal(extremeB[[1]]$value[2,2], mx, tolerance = 1e-10)
})
test_that("VARMA search works with predictions (mixture)", {
skip_on_cran()
res = search.varma(data = get.data(x[,1:7], endogenous = 3, newData = x[c(8,9,10),]),
combinations = get.combinations(sizes = c(1,2),
numTargets = 2,
innerGroups = list(c(1,2))),
maxParams = c(2,1,2,0,0,0),
maxHorizon = 3,
simUsePreviousEstim = FALSE,
metrics = get.search.metrics(c("sic"),c("rmse"),
horizons = c(1:3), simFixSize = 2),
items = get.search.items(all = TRUE, type1 = TRUE, mixture4 = TRUE),
modelChecks = get.search.modelchecks(estimation = FALSE, prediction = FALSE, predictionBound = 0))
sumRes <- summary(res, test = TRUE)
h = 2
coefs = c()
vars = c()
weights = c()
i <- 0
for (m in sumRes$results){
i <- i + 1
if (m$evalName != "rmse" || m$typeName != "model" || m$targetName != "V1")
next()
hh= h + m$value$prediction$startIndex - 1
coefs = append(coefs,m$value$prediction$means["V1",hh])
vars = append(vars, m$value$prediction$vars["V1",hh])
weights = append(weights, res$results[[i]]$value$weight)
}
mixture = res$results[which(sapply(res$results,
function(r) r$evalName == "rmse" && r$targetName == "V1" && r$typeName == "mixture"))]
# note that we need weighted mean, variance, etc. assuming normal distribution
me = weighted.mean(coefs, weights)
expect_equal(mixture[[1]]$value[2,1], me, tolerance = 1e-14)
# TODO : compare weighted variance, skewness, kurtosis assuming normality
# of course, its better to .Call the running statistics, test it, and use it here
#len = length(coefs)
# expect_equal(mixture[[1]]$value[2,5], len)
# test does not pass the second column for count (actual:48, res[2,2]: 46) in Debian
# note that mixture calculations are weight based
# same in CDF
})
test_that("estim.varma SplitSearch works (no subsetting)", {
skip_on_cran()
data = data = get.data(x[,1:7], endogenous = 3, newData = x[c(8,9,10),])
combinations = get.combinations(numTargets = 3, innerGroups = list(c(1), c(1,2), c(1,3)))
items = get.search.items(inclusion = TRUE
#, all = TRUE
, bestK = 4
, type1 = TRUE
, cdfs = c(0,0.3)
, mixture4 = TRUE
, extremeMultiplier = 2
)
metrics = get.search.metrics(c("sic", "aic"),
horizons = c(1:3), simFixSize = 2) # don't test with out-of-sample metrics. It seems we have different model with equal weights (the result change by repeating the call ?!)
options = get.search.options(FALSE,
#reportInterval = 1
)
combinations$sizes <- c(1, 2, 3)
whole = search.varma(data = data,
combinations = combinations,
maxParams = c(2,1,2,0,0,0),
maxHorizon = 3,
simUsePreviousEstim = FALSE,
items = items,
metrics = metrics,
options = options)
combinations$sizes <- list(c(1, 2), c(3))
combinations$stepsNumVariables <- c(NA, NA)
split = search.varma(data = data,
combinations = combinations,
maxParams = c(2,1,2,0,0,0),
maxHorizon = 3,
simUsePreviousEstim = FALSE,
items = items,
metrics = metrics,
options = options)
expect_equal(whole$counts, split$counts)
expect_equal(length(whole$results), length(split$results))
pastedList_w <- unlist(lapply(whole$results, function(x) paste(x[1:4], collapse = "")))
pastedList_s <- unlist(lapply(split$results, function(x) paste(x[1:4], collapse = "")))
sortedList_w <- whole$results[order(pastedList_w)]
sortedList_s <- split$results[order(pastedList_s)]
for (i in 1:length(sortedList_w)){
if (sortedList_s[[i]]$typeName == "mixture"){
expect_equal(sortedList_s[[i]]$value[,c(1:3,5,6)], sortedList_w[[i]]$value[,c(1:3,5,6)])
expect_equal(sortedList_s[[i]]$value[,c(4)], sortedList_w[[i]]$value[,c(4)], tolerance = 0.1)
}
else
expect_equal(sortedList_s[[i]]$value, sortedList_w[[i]]$value)
}
})
Try the ldt package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.