Nothing
tests/dse2tstgd2.R
In EvalEst: Dynamic Systems Estimation - Extensions
#######################################################################
# test functions for examples in the Brief User's Guide <<<<<<<<<<
#######################################################################
# The tests requiring padi and dsepadi are in dsepadi/tests.
if(!require("EvalEst")) stop("this test requires EvalEst.")
if(!require("setRNG"))stop("this test requires setRNG.")
Sys.info()
DSEversion()
print.values <- FALSE
fuzz.small <- 1e-14
fuzz.large <- 1e-8
max.error <- NA
all.ok <- TRUE
cat("Guide part 2 test 1 ... \n")
# Old Splus value
# {test.rng1 <- list(kind="default", normal.kind="default",
# seed=c(13,44,1,25,56,0,6,33,22,13,13,0) )
# test.rng2 <- list(kind="default", normal.kind="default",
# seed=c(13,43,7,57,62,3,30,29,24,54,47,2) )
# test.rng4 <- list(kind="default", normal.kind="default",
# seed=c(29,55,47,18,33,1,15,15,34,46,13,2) )
# test.rng3 <- list(kind="default", normal.kind="default",
# seed=c( 53,41,26,39,10,1,19,25,56,32,28,3) )
#As of R 1.7.1 a bug was recognized in Kinderman-Ramage and it was changed.
# These tests use the old version named "Buggy Kinderman-Ramage" to
# get old results for comparison.
test.rng1 <- test.rng2 <- test.rng3 <- test.rng4 <-
if (as.numeric(version$major)+0.1*as.numeric(version$minor) >= 1.71 )
list(kind="Wichmann-Hill", normal.kind="Buggy Kinderman-Ramage",
seed=c(979, 1479, 1542)) else
list(kind="Wichmann-Hill", normal.kind="Kinderman-Ramage",
seed=c(979, 1479, 1542))
# might consider also
#list(kind="Wichmann-Hill", normal.kind="user-supplied", seed=c(979, 1479, 1542))}
# R's Box-Muller was declared not reproducible.
#eg1.DSE.data <- t(matrix(scan(paste(DSE.HOME, "/data/eg1.dat", sep="")),
# 5, 364))[, 2:5]
#eg1.DSE.data <- TSdata(input = eg1.DSE.data[,1,drop = FALSE],
# output = eg1.DSE.data[, 2:4, drop = FALSE])
# above worked, but needs DSE.HOME which is depreciated. use
data("eg1.DSE.data", package="dse")
eg1.DSE.data <-tframed(eg1.DSE.data, list(start=c(1961,3), frequency=12))
seriesNamesInput(eg1.DSE.data) <- "R90"
seriesNamesOutput(eg1.DSE.data) <- c("M1","GDPl2", "CPI")
data("egJofF.1dec93.data", package="dse")
if(!exists("egJofF.1dec93.data"))warning("egJofF.1dec93.data does not exist")
error <- abs(3352.4721630925987483 -
sum(c(outputData(egJofF.1dec93.data),
inputData(egJofF.1dec93.data))))
ok <- fuzz.large > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok &ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
# Section 4 - Model Estimation
cat("Guide part 2 test 2 ... \n")
model.eg1.ls <- estVARXls(trimNA(eg1.DSE.data), warn=FALSE)
# opts <-options(warn=-1)
subsample.data <- tfwindow(eg1.DSE.data,start=c(1972,1),end=c(1992,12),warn=FALSE)
# options(opts)
# summary(model.eg1.ls)
# print(model.eg1.ls)
tfplot(model.eg1.ls)
tfplot(model.eg1.ls, start=c(1990,1))
z <- checkResiduals(model.eg1.ls, plot.=FALSE, pac=TRUE)
# Old Splus value c(4.67445135116577148, 3274.42578125, -2371.9997808950302)
check.value <- c(4.674448837156188, 3274.4226534894947, -2371.999780895034)
# above minor modification of the value in 1.8.0 beta
# else c(4.674448837156188, 3274.422653488969, -2371.999780895034) )
# using my old acf instead of bats version gives
# c(4.6744488371561879, 0.0, -2371.999780895033837)
tst <- c(sum(z$acf),sum(z$pacf),sum(z$cusum))
printTestValue(c(tst), digits=18)
error <- max(abs(check.value - tst ))
ok <- fuzz.large > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 3 ... \n")
# NB- non-stationary data. ar is not really valid
model.eg1.ar <- estVARXar(trimNA(eg1.DSE.data), warn=FALSE)
model.eg1.ss <- estSSfromVARX(trimNA(eg1.DSE.data), warn=FALSE)
# model.eg1.mle <- estMaxLik(trimNA(eg1.DSE.data),model.eg1.ar) # this may be slow
# Old Splus value c(6738.642280883833, 6921.352391513382)
check.value <- c(6738.562190977154, 6921.352391513382)#ts ar
# R bats ar c(6735.139112062216, 6921.352391513382)bats ar
# using my old ar:=ls gives c(6921.352391513380, 6921.352391513380)#ar:=ls
tst <- c(model.eg1.ar$estimates$like[1], model.eg1.ss$estimates$like[1])
printTestValue(c(tst), digits=18)
error <- max(abs(check.value - tst))
ok <- 10*fuzz.large > error
if ((Sys.info()[["sysname"]] == "Linux") && ! ok) {
warning("Using relaxed tolerance for Linux.")
ok <- 100*fuzz.large > error
}
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 4a... \n")
eg4.DSE.data<- egJofF.1dec93.data
outputData(eg4.DSE.data) <- outputData(eg4.DSE.data, series=c(1,2,6,7))
# following is optional
# tframe(outputData(eg4.DSE.data))<- tframe(outputData(egJofF.1dec93.data))
model.eg4.bb <- estBlackBox(trimNA(eg4.DSE.data), max.lag=3, verbose=FALSE)
tst <- model.eg4.bb$estimates$like[1]
printTestValue(c(tst), digits=18)
error <- abs(614.70500313590287078 - tst )
ok <- fuzz.large > error
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 4b... \n")
z <- informationTests(model.eg1.ar, model.eg1.ss, Print=FALSE, warn=FALSE)
# # Old Splus value check.value <- 231152.464267979725
# R check.value <- 231151.0300943982 # ts ar
# R check.value <- 230978.2532793634 bats ar
# R using my old ar=ls gives 233856.9237566061 #using ls for ar
# a small change in the accounting for degenerate subspaces in dse.2000.4 gives
# Splus value 225328.464267979754
check.value <- 225327.03009431256
if (print.values) printTestValue(sum(z[!is.na(z)]) )
error <- abs(check.value - sum(z[!is.na(z)]) )
ok <- 1000*fuzz.large > error #fuzz.large works in Solaris but not Linux
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
# Section 5 - Forecasting
cat("Guide part 2 test 5 pre-test a... \n")
eg4.DSE.model <- estVARXls(eg4.DSE.data)
longIn.data <- TSdata(
input= ts(rbind(inputData(eg4.DSE.data), matrix(.1,10,1)),
start=start(eg4.DSE.data),
frequency=frequency(eg4.DSE.data)),
output=outputData(eg4.DSE.data))
seriesNames(longIn.data) <- seriesNames(eg4.DSE.data)
z <- forecast(TSmodel(eg4.DSE.model), longIn.data)
zz <- forecast(TSmodel(eg4.DSE.model), longIn.data, compiled=FALSE)
error <- max(abs(z$pred - zz$pred))
ok <- fuzz.small > error
if (ok) cat("ok\n") else {
max.error <- if (is.na(max.error)) error else max(error, max.error)
cat("failed! error= ", error,"\n")
if(!testEqual(z$data$output, zz$data$output))
cat("output data comparison for forecast() and longIn.data failed.\n")
if(!testEqual(z$data$input, zz$data$input))
cat("input data comparison for forecast() and longIn.data failed.\n")
}
all.ok <- all.ok & ok
cat("Guide part 2 test 5 ... \n")
eg4.DSE.model <- estVARXls(eg4.DSE.data)
# Fame call disabled for testing: new.data <- freeze(eg4.DSE.data.names)
new.data <- TSdata(
input= ts(rbind(inputData(eg4.DSE.data), matrix(.1,10,1)),
start=start(eg4.DSE.data),
frequency=frequency(eg4.DSE.data)),
output=ts(rbind(outputData(eg4.DSE.data),matrix(.3,5,4)),
start=start(eg4.DSE.data),
frequency=frequency(eg4.DSE.data)))
seriesNames(new.data) <- seriesNames(eg4.DSE.data)
z <- l(TSmodel(eg4.DSE.model), trimNA(new.data))
# z <- l(TSmodel(eg4.DSE.model), trimNA(freeze(eg4.DSE.data.names)))
error <- max(abs(556.55870662521476788 -z$estimates$like[1]))
ok <- fuzz.large > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 6 ... \n")
zz <- forecast(TSmodel(eg4.DSE.model), new.data)
#tfplot(zz$pred, forecast(TSmodel(eg4.DSE.model), new.data, compiled=FALSE)$pred)
z <- forecast(TSmodel(eg4.DSE.model), trimNA(new.data),
conditioning.inputs=inputData(new.data))
tfplot(zz, start=c(1990,6))
error <- abs(4.7990339556773520258 - sum(forecasts(z)[[1]]))
ok <- fuzz.small > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
ok <- testEqual(zz,z) & ok
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 7 ... \n")
z <- forecast(eg4.DSE.model, conditioning.inputs.forecasts=matrix(.5,6,1))
# Fame call disabled for testing:
# z <- forecast(TSmodel(eg4.DSE.model), freeze(egJofF.1dec93.data.names),
# conditioning.inputs.forecasts=matrix(.5,6,1))
# summary(z)
# print(z)
tfplot(z)
tfplot(z, start=c(1990,1))
# forecasts(z)[[1]]
# tfwindow(forecasts(z)[[1]], start=c(1994,5))
ok <- all(start(eg4.DSE.model$data) == c(1974,2)) &
all(start(egJofF.1dec93.data) == c(1974,2))
error <- max(abs(c(5.9414711908521793404 - sum(forecasts(z)[[1]][1:6,]),
3.7410224783909828972 -
sum(tfwindow(forecasts(z)[[1]], start=c(1993,12), warn=FALSE)))))
ok <- ok & (fuzz.small > error)
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 8 ... \n")
z <- l(TSmodel(eg4.DSE.model), new.data)
tfplot(z)
z <- featherForecasts(TSmodel(eg4.DSE.model), new.data)
tfplot(z)
zz <-featherForecasts(TSmodel(eg4.DSE.model), new.data,
from.periods =c(20,50,60,70,80), horizon=150)
tfplot(zz)
error <- max(abs(c(54.838475604100473504 -
sum( forecasts(z)[[1]][10:46,]),
53.824873541066445171 -
sum(forecasts(zz)[[5]][80:116,]))))
ok <-fuzz.large > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 9 ... \n")
z <- horizonForecasts(TSmodel(eg4.DSE.model), new.data, horizons=c(1,3,6))
tfplot(z)
# error <- abs(653.329319170802592 - sum(z$horizonForecasts) )
error <- abs(653.329319170802592 - sum(forecasts(z)) )
ok <- fuzz.large > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 10... \n")
fc1 <- forecastCov(TSmodel(eg4.DSE.model), data=eg4.DSE.data)
tfplot(fc1)
tfplot(forecastCov(TSmodel(eg4.DSE.model), data=eg4.DSE.data, horizons= 1:4))
fc2 <- forecastCov(TSmodel(eg4.DSE.model),
data=eg4.DSE.data, zero=TRUE, trend=TRUE)
tfplot(fc2)
error <- max(abs(c(14.933660144821400806 - sum(fc1$forecastCov[[1]]),
14.933660144821400806 - sum(fc2$forecastCov[[1]]),
31.654672476928297442 - sum(fc2$forecastCov.zero),
18.324461923341953451 - sum(fc2$forecastCov.trend) )))
ok <- (fuzz.small * 1.5) > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
if (is.na(ok)) ok <- FALSE
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 11... \n")
mod1 <- ARMA(A=array(c(1,-.25,-.05), c(3,1,1)), B=array(1,c(1,1,1)))
mod2 <- ARMA(A=array(c(1,-.8, -.2 ), c(3,1,1)), B=array(1,c(1,1,1)))
mod3 <- ARMA(
A=array(c(
1.00,-0.06,0.15,-0.03,0.00,0.02,0.03,-0.02,0.00,-0.02,-0.03,-0.02,
0.00,-0.07,-0.05,0.12,1.00,0.20,-0.03,-0.11,0.00,-0.07,-0.03,0.08,
0.00,-0.40,-0.05,-0.66,0.00,0.00,0.17,-0.18,1.00,-0.11,-0.24,-0.09 )
,c(4,3,3)),
B=array(diag(1,3),c(1,3,3)))
e.ls.mod1 <- EstEval( mod1, replications=100,
rng=test.rng1,
simulation.args=list(sampleT=100, sd=1),
estimation="estVARXls", estimation.args=list(max.lag=2),
criterion="TSmodel", quiet=TRUE)
# e.ar.mod1 <- EstEval( mod1, replications=100,
# rng=test.rng1,
# simulation.args=list(sampleT=100, sd=1),
# estimation="estVARXar", estimation.args=list(max.lag=2, aic=FALSE),
# criterion="TSmodel", quiet=TRUE)
# tfplot(coef(e.ar.mod1))
# Old Splus value -0.29855874505752699744
check.value <- -0.3699580622977686
error <- abs( check.value - sum(coef(e.ls.mod1$result[[100]])))
ok <- fuzz.small > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 12... \n")
e.ls.mod2 <- EstEval( mod2, replications=100,
rng=test.rng2,
simulation.args=list(sampleT=100, sd=1),
estimation="estVARXls", estimation.args=list(max.lag=2),
criterion="TSmodel", quiet=TRUE)
old.par <- par(mfcol=c(2,1)) #set the number of plots on the plotics device
on.exit(par(old.par))
tfplot(coef(e.ls.mod1))
tfplot(coef(e.ls.mod2))
old.par <- c(old.par, par(mfcol=c(2,1)) )
tfplot(coef(e.ls.mod1), cum=FALSE, bounds=FALSE)
tfplot(coef(e.ls.mod2), cum=FALSE, bounds=FALSE)
distribution(coef(e.ls.mod1), bandwidth=.2)
distribution(coef(e.ls.mod2), bandwidth=.2)
# Old Splus value -1.0021490287427212706
check.value <- -1.0028944627996934
error <- abs(check.value - sum(coef(e.ls.mod2$result[[100]])))
ok <- fuzz.small > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 13... \n")
e.ls.mod1.roots <- roots(e.ls.mod1)
plot(e.ls.mod1.roots)
plot(e.ls.mod1.roots, complex.plane=FALSE)
plot(roots(e.ls.mod2), complex.plane=FALSE)
distribution(e.ls.mod1.roots, bandwidth=.2)
distribution(roots(e.ls.mod2), bandwidth=.1)
# Old Splus value 0.36159459310761993267
check.value <- 0.2119677206564640
# error <- Mod(0.36159459310761993267+0i - sum(e.ls.mod1.roots$result[[100]]))
error <- Mod(check.value - sum(e.ls.mod1.roots$result[[100]]))
ok <- fuzz.small > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 14... \n")
pc <- forecastCovEstimatorsWRTtrue(mod3,
rng=test.rng3,
estimation.methods=list(estVARXls=list(max.lag=6)),
est.replications=2, pred.replications=10, quiet=TRUE)
# the fuzz.small has to be relaxed here to accomodate differences in rnorm
# between Splus3.1 and Splus3.2 (the numbers are from Splus3.2)
# Old Splus value c(60.927013860328429473, 62.32729288591478678, 63.17808145947956433)
check.value <- c(54.164759056117504, 53.519297277839669, 59.341526159411558)
error <- max(abs(check.value -c(sum(pc$forecastCov[[1]]),
sum(pc$forecastCov.zero), sum(pc$forecastCov.trend) )))
ok <- fuzz.large > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 15... \n")
pc.rd <- forecastCovReductionsWRTtrue(mod3,
rng=test.rng4,
estimation.methods=list(estVARXls=list(max.lag=3)),
est.replications=2, pred.replications=10, quiet=TRUE)
# Old Splus value c(58.75543799264762157,60.451513998215133938, 64.089618782185240775)
check.value <- c(51.237201863944890, 53.519297277839669, 59.341526159411558)
error <- max(abs(check.value - c(sum(pc.rd$forecastCov[[1]]),
sum(pc.rd$forecastCov.zero), sum(pc.rd$forecastCov.trend))))
ok <- fuzz.large > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 16... \n")
z <-outOfSample.forecastCovEstimatorsWRTdata(trimNA(eg1.DSE.data),
estimation.sample=.5,
estimation.methods = list(
estVARXar=list(warn=FALSE),
estVARXls=list(warn=FALSE)),
trend=TRUE, zero=TRUE)
tfplot(z)
opts <- options(warn=-1)
zz <-outOfSample.forecastCovEstimatorsWRTdata(trimNA(eg1.DSE.data),
estimation.sample=.5,
estimation.methods = list(
estBlackBox4=list(max.lag=3, verbose=FALSE, warn=FALSE),
estVARXls=list(max.lag=3, warn=FALSE)),
trend=TRUE, zero=TRUE)
# zf<-horizonForecasts(zz$multi.model[[1]],zz$data, horizons=c(1,3,6))
zf<-horizonForecasts(TSmodel(zz, select=1),TSdata(zz), horizons=c(1,3,6))
options(opts)
zf<- zf$horizonForecasts[3,30,]
tfplot(z)
# Old Splus value
# c(6120.97621905043979, 175568.040899036743, 24.568074094041549,
# 1e-10*c(158049871127.845642, 3330592793.50789356,
# 1242727188.69001055, 1606263575.00784183))
check.value <- c(6120.97621905044, 175568.0408990367, 24.56807409404155,
1e-10*c(158034797997.4372, 3330592793.507894,
1242727188.690011, 1606263575.007842))
# using my old ls for ar instead of bats version gives
# c(6120.9762190509673, 175568.04089903546, 24.568074093999403,
# 1e-10*c(3330592793.507894, 3330592793.507894,
# 1242727188.690011, 1606263575.0078418)) # using ls for ar
if (print.values) printTestValue(c(zf, # 1e-5*
c(sum( z$forecastCov[[1]]), sum( z$forecastCov[[2]]),
sum(zz$forecastCov[[1]]), sum(zz$forecastCov[[2]]))) )
error <- max(abs(check.value -
c(zf, 1e-10*c(sum( z$forecastCov[[1]]), sum( z$forecastCov[[2]]),
sum(zz$forecastCov[[1]]), sum(zz$forecastCov[[2]])))))
ok <- fuzz.large > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("All Brief User Guide example tests part 2 completed")
if (all.ok) cat(" OK\n") else
cat(", some FAILED! max.error = ", max.error,"\n")
if (!all.ok) stop("Some tests FAILED")
Try the EvalEst package in your browser
Any scripts or data that you put into this service are public.
EvalEst documentation built on March 3, 2021, 1:14 a.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 functions for examples in the Brief User's Guide <<<<<<<<<<
#######################################################################
# The tests requiring padi and dsepadi are in dsepadi/tests.
if(!require("EvalEst")) stop("this test requires EvalEst.")
if(!require("setRNG"))stop("this test requires setRNG.")
Sys.info()
DSEversion()
print.values <- FALSE
fuzz.small <- 1e-14
fuzz.large <- 1e-8
max.error <- NA
all.ok <- TRUE
cat("Guide part 2 test 1 ... \n")
# Old Splus value
# {test.rng1 <- list(kind="default", normal.kind="default",
# seed=c(13,44,1,25,56,0,6,33,22,13,13,0) )
# test.rng2 <- list(kind="default", normal.kind="default",
# seed=c(13,43,7,57,62,3,30,29,24,54,47,2) )
# test.rng4 <- list(kind="default", normal.kind="default",
# seed=c(29,55,47,18,33,1,15,15,34,46,13,2) )
# test.rng3 <- list(kind="default", normal.kind="default",
# seed=c( 53,41,26,39,10,1,19,25,56,32,28,3) )
#As of R 1.7.1 a bug was recognized in Kinderman-Ramage and it was changed.
# These tests use the old version named "Buggy Kinderman-Ramage" to
# get old results for comparison.
test.rng1 <- test.rng2 <- test.rng3 <- test.rng4 <-
if (as.numeric(version$major)+0.1*as.numeric(version$minor) >= 1.71 )
list(kind="Wichmann-Hill", normal.kind="Buggy Kinderman-Ramage",
seed=c(979, 1479, 1542)) else
list(kind="Wichmann-Hill", normal.kind="Kinderman-Ramage",
seed=c(979, 1479, 1542))
# might consider also
#list(kind="Wichmann-Hill", normal.kind="user-supplied", seed=c(979, 1479, 1542))}
# R's Box-Muller was declared not reproducible.
#eg1.DSE.data <- t(matrix(scan(paste(DSE.HOME, "/data/eg1.dat", sep="")),
# 5, 364))[, 2:5]
#eg1.DSE.data <- TSdata(input = eg1.DSE.data[,1,drop = FALSE],
# output = eg1.DSE.data[, 2:4, drop = FALSE])
# above worked, but needs DSE.HOME which is depreciated. use
data("eg1.DSE.data", package="dse")
eg1.DSE.data <-tframed(eg1.DSE.data, list(start=c(1961,3), frequency=12))
seriesNamesInput(eg1.DSE.data) <- "R90"
seriesNamesOutput(eg1.DSE.data) <- c("M1","GDPl2", "CPI")
data("egJofF.1dec93.data", package="dse")
if(!exists("egJofF.1dec93.data"))warning("egJofF.1dec93.data does not exist")
error <- abs(3352.4721630925987483 -
sum(c(outputData(egJofF.1dec93.data),
inputData(egJofF.1dec93.data))))
ok <- fuzz.large > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok &ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
# Section 4 - Model Estimation
cat("Guide part 2 test 2 ... \n")
model.eg1.ls <- estVARXls(trimNA(eg1.DSE.data), warn=FALSE)
# opts <-options(warn=-1)
subsample.data <- tfwindow(eg1.DSE.data,start=c(1972,1),end=c(1992,12),warn=FALSE)
# options(opts)
# summary(model.eg1.ls)
# print(model.eg1.ls)
tfplot(model.eg1.ls)
tfplot(model.eg1.ls, start=c(1990,1))
z <- checkResiduals(model.eg1.ls, plot.=FALSE, pac=TRUE)
# Old Splus value c(4.67445135116577148, 3274.42578125, -2371.9997808950302)
check.value <- c(4.674448837156188, 3274.4226534894947, -2371.999780895034)
# above minor modification of the value in 1.8.0 beta
# else c(4.674448837156188, 3274.422653488969, -2371.999780895034) )
# using my old acf instead of bats version gives
# c(4.6744488371561879, 0.0, -2371.999780895033837)
tst <- c(sum(z$acf),sum(z$pacf),sum(z$cusum))
printTestValue(c(tst), digits=18)
error <- max(abs(check.value - tst ))
ok <- fuzz.large > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 3 ... \n")
# NB- non-stationary data. ar is not really valid
model.eg1.ar <- estVARXar(trimNA(eg1.DSE.data), warn=FALSE)
model.eg1.ss <- estSSfromVARX(trimNA(eg1.DSE.data), warn=FALSE)
# model.eg1.mle <- estMaxLik(trimNA(eg1.DSE.data),model.eg1.ar) # this may be slow
# Old Splus value c(6738.642280883833, 6921.352391513382)
check.value <- c(6738.562190977154, 6921.352391513382)#ts ar
# R bats ar c(6735.139112062216, 6921.352391513382)bats ar
# using my old ar:=ls gives c(6921.352391513380, 6921.352391513380)#ar:=ls
tst <- c(model.eg1.ar$estimates$like[1], model.eg1.ss$estimates$like[1])
printTestValue(c(tst), digits=18)
error <- max(abs(check.value - tst))
ok <- 10*fuzz.large > error
if ((Sys.info()[["sysname"]] == "Linux") && ! ok) {
warning("Using relaxed tolerance for Linux.")
ok <- 100*fuzz.large > error
}
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 4a... \n")
eg4.DSE.data<- egJofF.1dec93.data
outputData(eg4.DSE.data) <- outputData(eg4.DSE.data, series=c(1,2,6,7))
# following is optional
# tframe(outputData(eg4.DSE.data))<- tframe(outputData(egJofF.1dec93.data))
model.eg4.bb <- estBlackBox(trimNA(eg4.DSE.data), max.lag=3, verbose=FALSE)
tst <- model.eg4.bb$estimates$like[1]
printTestValue(c(tst), digits=18)
error <- abs(614.70500313590287078 - tst )
ok <- fuzz.large > error
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 4b... \n")
z <- informationTests(model.eg1.ar, model.eg1.ss, Print=FALSE, warn=FALSE)
# # Old Splus value check.value <- 231152.464267979725
# R check.value <- 231151.0300943982 # ts ar
# R check.value <- 230978.2532793634 bats ar
# R using my old ar=ls gives 233856.9237566061 #using ls for ar
# a small change in the accounting for degenerate subspaces in dse.2000.4 gives
# Splus value 225328.464267979754
check.value <- 225327.03009431256
if (print.values) printTestValue(sum(z[!is.na(z)]) )
error <- abs(check.value - sum(z[!is.na(z)]) )
ok <- 1000*fuzz.large > error #fuzz.large works in Solaris but not Linux
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
# Section 5 - Forecasting
cat("Guide part 2 test 5 pre-test a... \n")
eg4.DSE.model <- estVARXls(eg4.DSE.data)
longIn.data <- TSdata(
input= ts(rbind(inputData(eg4.DSE.data), matrix(.1,10,1)),
start=start(eg4.DSE.data),
frequency=frequency(eg4.DSE.data)),
output=outputData(eg4.DSE.data))
seriesNames(longIn.data) <- seriesNames(eg4.DSE.data)
z <- forecast(TSmodel(eg4.DSE.model), longIn.data)
zz <- forecast(TSmodel(eg4.DSE.model), longIn.data, compiled=FALSE)
error <- max(abs(z$pred - zz$pred))
ok <- fuzz.small > error
if (ok) cat("ok\n") else {
max.error <- if (is.na(max.error)) error else max(error, max.error)
cat("failed! error= ", error,"\n")
if(!testEqual(z$data$output, zz$data$output))
cat("output data comparison for forecast() and longIn.data failed.\n")
if(!testEqual(z$data$input, zz$data$input))
cat("input data comparison for forecast() and longIn.data failed.\n")
}
all.ok <- all.ok & ok
cat("Guide part 2 test 5 ... \n")
eg4.DSE.model <- estVARXls(eg4.DSE.data)
# Fame call disabled for testing: new.data <- freeze(eg4.DSE.data.names)
new.data <- TSdata(
input= ts(rbind(inputData(eg4.DSE.data), matrix(.1,10,1)),
start=start(eg4.DSE.data),
frequency=frequency(eg4.DSE.data)),
output=ts(rbind(outputData(eg4.DSE.data),matrix(.3,5,4)),
start=start(eg4.DSE.data),
frequency=frequency(eg4.DSE.data)))
seriesNames(new.data) <- seriesNames(eg4.DSE.data)
z <- l(TSmodel(eg4.DSE.model), trimNA(new.data))
# z <- l(TSmodel(eg4.DSE.model), trimNA(freeze(eg4.DSE.data.names)))
error <- max(abs(556.55870662521476788 -z$estimates$like[1]))
ok <- fuzz.large > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 6 ... \n")
zz <- forecast(TSmodel(eg4.DSE.model), new.data)
#tfplot(zz$pred, forecast(TSmodel(eg4.DSE.model), new.data, compiled=FALSE)$pred)
z <- forecast(TSmodel(eg4.DSE.model), trimNA(new.data),
conditioning.inputs=inputData(new.data))
tfplot(zz, start=c(1990,6))
error <- abs(4.7990339556773520258 - sum(forecasts(z)[[1]]))
ok <- fuzz.small > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
ok <- testEqual(zz,z) & ok
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 7 ... \n")
z <- forecast(eg4.DSE.model, conditioning.inputs.forecasts=matrix(.5,6,1))
# Fame call disabled for testing:
# z <- forecast(TSmodel(eg4.DSE.model), freeze(egJofF.1dec93.data.names),
# conditioning.inputs.forecasts=matrix(.5,6,1))
# summary(z)
# print(z)
tfplot(z)
tfplot(z, start=c(1990,1))
# forecasts(z)[[1]]
# tfwindow(forecasts(z)[[1]], start=c(1994,5))
ok <- all(start(eg4.DSE.model$data) == c(1974,2)) &
all(start(egJofF.1dec93.data) == c(1974,2))
error <- max(abs(c(5.9414711908521793404 - sum(forecasts(z)[[1]][1:6,]),
3.7410224783909828972 -
sum(tfwindow(forecasts(z)[[1]], start=c(1993,12), warn=FALSE)))))
ok <- ok & (fuzz.small > error)
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 8 ... \n")
z <- l(TSmodel(eg4.DSE.model), new.data)
tfplot(z)
z <- featherForecasts(TSmodel(eg4.DSE.model), new.data)
tfplot(z)
zz <-featherForecasts(TSmodel(eg4.DSE.model), new.data,
from.periods =c(20,50,60,70,80), horizon=150)
tfplot(zz)
error <- max(abs(c(54.838475604100473504 -
sum( forecasts(z)[[1]][10:46,]),
53.824873541066445171 -
sum(forecasts(zz)[[5]][80:116,]))))
ok <-fuzz.large > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 9 ... \n")
z <- horizonForecasts(TSmodel(eg4.DSE.model), new.data, horizons=c(1,3,6))
tfplot(z)
# error <- abs(653.329319170802592 - sum(z$horizonForecasts) )
error <- abs(653.329319170802592 - sum(forecasts(z)) )
ok <- fuzz.large > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 10... \n")
fc1 <- forecastCov(TSmodel(eg4.DSE.model), data=eg4.DSE.data)
tfplot(fc1)
tfplot(forecastCov(TSmodel(eg4.DSE.model), data=eg4.DSE.data, horizons= 1:4))
fc2 <- forecastCov(TSmodel(eg4.DSE.model),
data=eg4.DSE.data, zero=TRUE, trend=TRUE)
tfplot(fc2)
error <- max(abs(c(14.933660144821400806 - sum(fc1$forecastCov[[1]]),
14.933660144821400806 - sum(fc2$forecastCov[[1]]),
31.654672476928297442 - sum(fc2$forecastCov.zero),
18.324461923341953451 - sum(fc2$forecastCov.trend) )))
ok <- (fuzz.small * 1.5) > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
if (is.na(ok)) ok <- FALSE
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 11... \n")
mod1 <- ARMA(A=array(c(1,-.25,-.05), c(3,1,1)), B=array(1,c(1,1,1)))
mod2 <- ARMA(A=array(c(1,-.8, -.2 ), c(3,1,1)), B=array(1,c(1,1,1)))
mod3 <- ARMA(
A=array(c(
1.00,-0.06,0.15,-0.03,0.00,0.02,0.03,-0.02,0.00,-0.02,-0.03,-0.02,
0.00,-0.07,-0.05,0.12,1.00,0.20,-0.03,-0.11,0.00,-0.07,-0.03,0.08,
0.00,-0.40,-0.05,-0.66,0.00,0.00,0.17,-0.18,1.00,-0.11,-0.24,-0.09 )
,c(4,3,3)),
B=array(diag(1,3),c(1,3,3)))
e.ls.mod1 <- EstEval( mod1, replications=100,
rng=test.rng1,
simulation.args=list(sampleT=100, sd=1),
estimation="estVARXls", estimation.args=list(max.lag=2),
criterion="TSmodel", quiet=TRUE)
# e.ar.mod1 <- EstEval( mod1, replications=100,
# rng=test.rng1,
# simulation.args=list(sampleT=100, sd=1),
# estimation="estVARXar", estimation.args=list(max.lag=2, aic=FALSE),
# criterion="TSmodel", quiet=TRUE)
# tfplot(coef(e.ar.mod1))
# Old Splus value -0.29855874505752699744
check.value <- -0.3699580622977686
error <- abs( check.value - sum(coef(e.ls.mod1$result[[100]])))
ok <- fuzz.small > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 12... \n")
e.ls.mod2 <- EstEval( mod2, replications=100,
rng=test.rng2,
simulation.args=list(sampleT=100, sd=1),
estimation="estVARXls", estimation.args=list(max.lag=2),
criterion="TSmodel", quiet=TRUE)
old.par <- par(mfcol=c(2,1)) #set the number of plots on the plotics device
on.exit(par(old.par))
tfplot(coef(e.ls.mod1))
tfplot(coef(e.ls.mod2))
old.par <- c(old.par, par(mfcol=c(2,1)) )
tfplot(coef(e.ls.mod1), cum=FALSE, bounds=FALSE)
tfplot(coef(e.ls.mod2), cum=FALSE, bounds=FALSE)
distribution(coef(e.ls.mod1), bandwidth=.2)
distribution(coef(e.ls.mod2), bandwidth=.2)
# Old Splus value -1.0021490287427212706
check.value <- -1.0028944627996934
error <- abs(check.value - sum(coef(e.ls.mod2$result[[100]])))
ok <- fuzz.small > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 13... \n")
e.ls.mod1.roots <- roots(e.ls.mod1)
plot(e.ls.mod1.roots)
plot(e.ls.mod1.roots, complex.plane=FALSE)
plot(roots(e.ls.mod2), complex.plane=FALSE)
distribution(e.ls.mod1.roots, bandwidth=.2)
distribution(roots(e.ls.mod2), bandwidth=.1)
# Old Splus value 0.36159459310761993267
check.value <- 0.2119677206564640
# error <- Mod(0.36159459310761993267+0i - sum(e.ls.mod1.roots$result[[100]]))
error <- Mod(check.value - sum(e.ls.mod1.roots$result[[100]]))
ok <- fuzz.small > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 14... \n")
pc <- forecastCovEstimatorsWRTtrue(mod3,
rng=test.rng3,
estimation.methods=list(estVARXls=list(max.lag=6)),
est.replications=2, pred.replications=10, quiet=TRUE)
# the fuzz.small has to be relaxed here to accomodate differences in rnorm
# between Splus3.1 and Splus3.2 (the numbers are from Splus3.2)
# Old Splus value c(60.927013860328429473, 62.32729288591478678, 63.17808145947956433)
check.value <- c(54.164759056117504, 53.519297277839669, 59.341526159411558)
error <- max(abs(check.value -c(sum(pc$forecastCov[[1]]),
sum(pc$forecastCov.zero), sum(pc$forecastCov.trend) )))
ok <- fuzz.large > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 15... \n")
pc.rd <- forecastCovReductionsWRTtrue(mod3,
rng=test.rng4,
estimation.methods=list(estVARXls=list(max.lag=3)),
est.replications=2, pred.replications=10, quiet=TRUE)
# Old Splus value c(58.75543799264762157,60.451513998215133938, 64.089618782185240775)
check.value <- c(51.237201863944890, 53.519297277839669, 59.341526159411558)
error <- max(abs(check.value - c(sum(pc.rd$forecastCov[[1]]),
sum(pc.rd$forecastCov.zero), sum(pc.rd$forecastCov.trend))))
ok <- fuzz.large > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("Guide part 2 test 16... \n")
z <-outOfSample.forecastCovEstimatorsWRTdata(trimNA(eg1.DSE.data),
estimation.sample=.5,
estimation.methods = list(
estVARXar=list(warn=FALSE),
estVARXls=list(warn=FALSE)),
trend=TRUE, zero=TRUE)
tfplot(z)
opts <- options(warn=-1)
zz <-outOfSample.forecastCovEstimatorsWRTdata(trimNA(eg1.DSE.data),
estimation.sample=.5,
estimation.methods = list(
estBlackBox4=list(max.lag=3, verbose=FALSE, warn=FALSE),
estVARXls=list(max.lag=3, warn=FALSE)),
trend=TRUE, zero=TRUE)
# zf<-horizonForecasts(zz$multi.model[[1]],zz$data, horizons=c(1,3,6))
zf<-horizonForecasts(TSmodel(zz, select=1),TSdata(zz), horizons=c(1,3,6))
options(opts)
zf<- zf$horizonForecasts[3,30,]
tfplot(z)
# Old Splus value
# c(6120.97621905043979, 175568.040899036743, 24.568074094041549,
# 1e-10*c(158049871127.845642, 3330592793.50789356,
# 1242727188.69001055, 1606263575.00784183))
check.value <- c(6120.97621905044, 175568.0408990367, 24.56807409404155,
1e-10*c(158034797997.4372, 3330592793.507894,
1242727188.690011, 1606263575.007842))
# using my old ls for ar instead of bats version gives
# c(6120.9762190509673, 175568.04089903546, 24.568074093999403,
# 1e-10*c(3330592793.507894, 3330592793.507894,
# 1242727188.690011, 1606263575.0078418)) # using ls for ar
if (print.values) printTestValue(c(zf, # 1e-5*
c(sum( z$forecastCov[[1]]), sum( z$forecastCov[[2]]),
sum(zz$forecastCov[[1]]), sum(zz$forecastCov[[2]]))) )
error <- max(abs(check.value -
c(zf, 1e-10*c(sum( z$forecastCov[[1]]), sum( z$forecastCov[[2]]),
sum(zz$forecastCov[[1]]), sum(zz$forecastCov[[2]])))))
ok <- fuzz.large > error
if (!ok) {if (is.na(max.error)) max.error <- error
else max.error <- max(error, max.error)}
all.ok <- all.ok & ok
{if (ok) cat("ok\n") else cat("failed! error= ", error,"\n") }
cat("All Brief User Guide example tests part 2 completed")
if (all.ok) cat(" OK\n") else
cat(", some FAILED! max.error = ", max.error,"\n")
if (!all.ok) stop("Some tests FAILED")
Try the EvalEst 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.