mingstat / zim Goto Github PK

Hi, thanks for writing and making your package publicly available!

If I understand correct there is no straightforward way to forecast with ZIM right now. I'd like to make a request for a forecast (or predict) method.

Since the forecast package already includes a broad range of univariate forecasting models with a consistent API, it would be nice to have something that is easy to use in a similar fashion, e.g.:

library("forecast")
library("ZIM")

data(syph)
count <- syph$a33
count.lag1 <- bshift(count > 0)
trend <- 1:length(count) / 1000

mdl1 <- auto.arima(count)
forecast(mdl1, h = 3)
# Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
# 210       2.836695 -1.014251 6.687642 -3.052820 8.726211
# 211       2.636530 -1.243722 6.516783 -3.297804 8.570865
# 212       2.618010 -1.265695 6.501715 -3.321605 8.557624

# with covariates
mdl2 <- auto.arima(head(count, -3), xreg = head(trend, -3))
forecast(mdl2, h = 3, xreg = tail(trend, 3))
# Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
# 207       2.148899 -1.655740 5.953538 -3.669794 7.967592
# 208       2.135985 -1.668653 5.940624 -3.682707 7.954678
# 209       2.123072 -1.681567 5.927711 -3.695621 7.941765

# ZIM
mdl3 <- zim(count ~ count.lag1 + trend | trend)
# e.g.:
forecast(mdl1, h = 6, xreg = ..., zreg = ...)
# something something here

It appears to be something wrong with an if defined inside "zim" method.
Next appears the code which gives the error...

count = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 650.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 300.0, 0.0, 0.0, 0.0, 0.0, 2100.0, 1600.0, 400.0, 750.0, 0.0, 0.0, 0.0, 1400.0, 620.0, 0.0, 600.0, 600.0, 0.0, 0.0, 2300.0, 0.0, 1950.0, 0.0, 750.0, 0.0, 0.0, 500.0, 300.0, 0.0, 0.0, 0.0, 0.0, 0.0, 350.0, 0.0, 350.0, 0.0, 2600.0, 0.0, 0.0, 2000.0, 0.0, 750.0, 5500.0, 900.0, 0.0, 0.0, 0.0, 350.0, 2200.0, 200.0, 0.0, 0.0, 0.0, 480.0, 0.0, 0.0, 0.0, 350.0, 0.0, 0.0, 340.0, 200.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1700.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 440.0, 300.0, 0.0, 0.0, 140.0, 200.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1000.0, 350.0, 190.0, 0.0, 400.0, 0.0, 0.0, 0.0, 900.0, 750.0, 0.0, 850.0, 0.0, 0.0, 0.0, 1000.0, 0.0, 800.0, 300.0, 0.0, 0.0, 590.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1680.0, 1000.0, 350.0, 0.0, 0.0, 0.0, 0.0, 475.9, 500.0, 2500.0, 200.0, 200.0, 0.0, 0.0, 460.0, 0.0, 2700.0, 0.0, 1800.0, 0.0, 0.0, 1600.0, 450.0, 0.0, 0.0, 500.0, 0.0, 0.0, 1100.0, 0.0, 0.0, 250.0, 0.0, 0.0, 0.0, 600.0, 2800.0, 0.0, 1160.0, 0.0, 0.0, 0.0, 3000.0, 750.0, 0.0, 350.0, 1900.0, 0.0, 0.0, 0.0, 0.0, 2200.0, 0.0, 0.0, 0.0, 0.0, 600.0, 0.0, 0.0, 0.0, 900.0, 0.0, 0.0, 600.0, 0.0, 0.0, 500.0, 0.0, 0.0, 0.0, 0.0, 0.0, 880.0, 0.0, 0.0, 0.0, 0.0, 800.0, 0.0, 0.0, 3290.0, 0.0, 0.0, 0.0, 1150.0, 0.0, 3000.0, 0.0, 300.0, 0.0, 0.0, 2100.0, 0.0, 600.0, 0.0, 1300.0, 0.0, 0.0, 0.0, 90.0, 150.0, 0.0, 1000.0, 0.0, 0.0, 1400.0, 0.0, 0.0, 0.0, 2800.0, 0.0, 0.0, 0.0, 1100.0, 0.0, 1150.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 600.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2900.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2600.0, 0.0, 0.0, 0.0, 0.0, 1050.0, 1600.0, 0.0, 0.0, 0.0, 0.0, 0.0, 600.0, 0.0, 0.0, 1300.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1650.0, 0.0, 0.0, 0.0, 850.0, 0.0, 1350.0, 0.0, 0.0, 0.0, 0.0, 900.0, 0.0, 0.0, 0.0, 1150.0, 0.0, 0.0, 0.0, 0.0, 900.0, 490.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1200.0, 0.0, 0.0, 0.0, 0.0, 1700.0, 2550.0, 0.0)
count = floor(count)
count.lag1 <- bshift(count > 0)
trend <- 1:length(count) / 1000
m1 = zim(count ~ count.lag1 + trend | trend)
Error in if (abs(loglik.new - loglik.old)/(abs(loglik.new) + control$epsilon) > :
valor ausente donde TRUE/FALSE es necesario
Translation: (missing value where TRUE/FALSE is needed)

If I use the negative binomial distribution all goes well.