Update fourier and fourierf to handle msts objects.

Update Arima and auto.arima to handle msts objects using the ts.frequency attribute.

Update msts so ts.frequency is not over-ridden by the existing frequency when argument is a ts object.

Add functions to compute trading day variables, Easter variable (as in Bell and Hillmer, JASA 1983), etc.

Hello Mr. Hyndman,

one of the initializations in your ets method seems weird:
https://github.com/robjhyndman/forecast/blob/master/R/ets.R#L622

In case n < 3m .. lets assume an additive season:
The mean of the first m periods is subtracted from the whole time series. Then you subtract these differences from the time series again .. which results in the mean of course. So now it's a straight line.
A few lines later a linear model (lsfit) is fitted where the intercept is of course .. the mean and the slope is nearly zero. Don't see the point of it. And by the way: why are the seasonal indices reverted but only from [2:m]?

Best regards.

For use in accuracy(), forecasterrors(), and elsewhere.

getResponse.default <- function(object){object$x}

This means other packages can use accuracy() even if x is not an element of the list. All they would need to do is define a getResponse method for the relevant class.

The command

simulate(auto.arima(as.vector(USAccDeaths)), 36)

yields Error in is.na(x) : 'x' is missing

but the example in simulate.Arima works perfectly fine (not using auto.arima)
Version 4.01

If all or nearly all observations in a time series are zero, Croston will fail. Example:

## Works fine
> croston(c(1,0,0,0,2,0,0), h=5, alpha=0.1)
   Point Forecast
8       0.8461538
9       0.8461538
10      0.8461538
11      0.8461538
12      0.8461538

## Failure with a single non-zero
> croston(c(0,0,0,0,2,0,0), h=5, alpha=0.1)
Error in matrix(x$mean, ncol = 1) : 'data' must be of a vector type
In addition: Warning message:
In croston(c(0, 0, 0, 0, 2, 0, 0), h = 5, alpha = 0.1) :
  Coercing LHS to a list

## Failure with all zeros
> croston(c(0,0,0,0,0,0,0), h=5, alpha=0.1)
Error in matrix(x$mean, ncol = 1) : 'data' must be of a vector type
In addition: Warning message:
In croston(c(0, 0, 0, 0, 0, 0, 0), h = 5, alpha = 0.1) :
  Coercing LHS to a list

This is with forecast 4.0.4 (latest on CRAN as of this post). I quickly stepped through the code and it's basically an edge case when doing the smoothing on the series of non-zeroes.

Thanks!

• Nicolas

When using xreg and newxreg, auto.arima can be very slow even with stepwise. Try using auto.arima with Canada[,1] from vars package and Canada[, 2:3] as xreg. Can this be speeded up?

See http://stats.stackexchange.com/q/76761/159

The following code is a reproducible example:

tser <- ts(c(0,0,0,0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,64,0,0),start=c(2010,1),freq=12)
mod.aach = auto.arima(tser, seasonal.test="ch")

The following error shows up:

Error in solve.default(t(A) %*% Omfhat %*% A, tol = 1e-25) :  system is computationally singular: reciprocal condition number = 6.50014e-31

If I try to use seasonal.test="ocsb" the model is successfully computed.

Is it an expected behaviour or something that can be fixed?

Currently not allowed, but would be possible by turning off Box-Cox transformations. In the case of zeros, a Box-Cox transform with lambda>0 would also work. Need to make it more flexible.

I'm not sure why these are so slow.

See http://stats.stackexchange.com/q/72244/159

library(tseries)
library(forecast)
dwa <- get.hist.quote(instrument="DWA", start="2010-01-01", end="2013-10-31")
logreturns <- diff(log(dwa$Close))*100
fit <- auto.arima(logreturns, trace=TRUE)

logreturns_nots = unclass(logreturns) ## converst time series into a vector.
fit <- auto.arima(logreturns_nots, trace=TRUE)

Either capture with an error message, or convert to ts class. This is ugly.

s <- ts(read.csv("Sample.csv",header=FALSE))
model <- tbats(s)
forecast(model, h=12)

Inf forecasts.

Hello sir,
I'm an undergraduate student using demography package to study Mauritian population. I'm having some problem using pop.sim function. How must the mid year population be entered?
Can you give me an example please?
Thank you very much in advance

Allow arfima() to take a model argument as with Arima and ets

Hello,

I want to get "Error Measures" (ME, RMSE,...) from output of summary(forecast( ... )) into a variable , so I can process them. it doesn't work ? There is no name for them in the forecast object !

R version 3.0.2 (2013-09-25) -- "Frisbee Sailing"
Copyright (C) 2013 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

library(forecast)
This is forecast 4.8

fit <- auto.arima(WWWusage)
a=forecast(fit,h=2)
a
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
101 219.1572 215.1402 223.1743 213.0137 225.3008
102 219.0323 209.4558 228.6087 204.3864 233.6782
summary(a)

Forecast method: ARIMA(1,1,1) with drift

Model Information:
Series: WWWusage
ARIMA(1,1,1) with drift

Coefficients:
ar1 ma1 drift
0.6344 0.5297 1.1205
s.e. 0.0871 0.0897 1.2925

sigma^2 estimated as 9.825: log likelihood=-251.23
AIC=510.45 AICc=510.88 BIC=520.83

Error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.04803303 3.10304 2.395923 0.073378 1.914083 0.5294563 -0.007904884

Forecasts:
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
101 219.1572 215.1402 223.1743 213.0137 225.3008
102 219.0323 209.4558 228.6087 204.3864 233.6782

names(a)
[1] "method" "model" "level" "mean" "lower" "upper" "x" "xname"
[9] "fitted" "residuals"

train <- window(a10, start = c(1991, 7), end = c(2007, 12))
test <- window(a10, start = c(2008, 1))
fcasts <- forecast(train, h = 6)
accuracy(fcasts, test, test = 1:2) # Warning (not expected)
accuracy(fcasts, test, test = 1:6)

Should return the same results as test = "all", but ACF1 is missing

accuracy(fcasts, test, test = “all")
accuracy(fcasts, test, test = 1:10) # Needs a proper warning indicating test exceeds indices
accuracy(fcasts, test, test = 2) # Error. cannot test a single element?

I would like an arar() function, or better still an ararma function().

See Newton-Parzen for details.

There is an arar() function already in the itsmr package.

See forecast.ets() for analogous functionality already available for forecast.ets().

Hi Rob,

It works fine when lambda = 0 and PI = TRUE, but not for PI = FALSE.

set.seed(1234)
data <- rpois(50, 20)
fit <- ets(data, lambda = 0)
pred <- forecast(fit, h = 4, lambda = 0, PI = FALSE)
## Error in exp(x) : non-numeric argument to mathematical function

ets() optimizes the initial states, but traditional exponential smoothing does not.
HoltWinters() does not optimize the initial states, but it uses a non-standard heuristic initialization.
For teaching purposes, it would be good to have functions that did not optimize initial states but used standard heuristic initialization.

I envisage functions ses(), holt() and hw(), each with an argument "initial" taking values "optimal" or "heuristic". If initial=optimal, the function will just call ets(). Otherwise it will do something like HoltWinters() but with more standard initialization methods.

If the time series has heteroscedasticity (http://en.wikipedia.org/wiki/Heteroscedasticity) should it be fixed before running auto.arima or will auto.arima detect and fix automatically?

I'm trying the new tsclean function with a univariate time series, with daily frequency, and I get the following error when calling it as tsclean(ts(data,frequency=7)):

Error in stl(xx, s.window = "periodic", robust = TRUE) :
only univariate series are allowed

Nevertheless, when I call the function as tsclean(ts(data)), it works (i.e., no error appears).

Cheers,

Mila

I realize this might not be a 'forecast' issue so much as a 'caret' issue, but maybe you can forward onto the maintainer for that package. Forecast depends on caret which depends on Hadley Wickam's 'reshape'. However I believe reshape is no longer updated and the preferred package is the much faster and more flexible reshape2. This is an issue for some because reshape masks a lot of the key functions for reshape2 which in turn can break a lot of existing code (at least it did for me). Again, I know this not really a 'forecast' problem but I thought I would bring it to your attention since I believe alot of people rely on both forecast and reshape2. Also, since reshape2 is notably faster than reshape, changing the dependence might bring about some performance upgrades. Thanks for all the great work.

Addendum: I contacted the caret package maintainer (Max) and he said this issue will be addressed in the next update (around late December).

This is for when users don't realise they already have the forecasts. e.g., with the dshw() function.

forecast.forecast <- function(object, ...)
{
return(object)
}

Could make it fancier by allowing different h and different level arguments.

Currently (v4.9), forecast.Arima(...,bootstrap=FALSE) still assumes normally distributed residuals. This should be changed to compute percentiles from future sample paths based on bootstrapped residuals.

Want feed-forward neural network with one hidden layer.
User to specify embedding dimension and size.
Need forecast method as well.

Allow a specific model to be specified in both bats and tbats.
Also allow a fitted model to be reused, as in ets and Arima.

I have two machines with the exactly same version of R. One can run the following code, the other gives an error.

vec <- c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1045.0, 872.0, 824.0, 1148.0, 1222.0, 940.0, 727.0, 770.0, 879.0, 960.0)
tser <- ts(vec,start=c(2010,1),freq=12)
auto.arima(tser,seasonal.test="ch")

This is the error I got in one machine:

Error in solve.default(tmp, tol = 1e-25) : 
  Lapack routine dgesv: system is exactly singular: U[11,11] = 0

INFO MACHINE N.1 (NO ERROR)

OS INFO

Distributor ID: Ubuntu
Description: Ubuntu 12.04.2 LTS
Release: 12.04
Codename: precise

R INFO

platform x86_64-pc-linux-gnu
arch x86_64
os linux-gnu
system x86_64, linux-gnu
status
major 3
minor 0.1
year 2013
month 05
day 16
svn rev 62743
language R
version.string R version 3.0.1 (2013-05-16)
nickname Good Sport

INFO MACHINE N.2 (WITH ERROR)

OS INFO

Distributor ID: Ubuntu
Description: Ubuntu 12.04.3 LTS
Release: 12.04
Codename: precise

R INFO

platform x86_64-pc-linux-gnu
arch x86_64
os linux-gnu
system x86_64, linux-gnu
status
major 3
minor 0.1
year 2013
month 05
day 16
svn rev 62743
language R
version.string R version 3.0.1 (2013-05-16)
nickname Good Sport

Resuming, the only difference is the third level version of Ubuntu.

Another thing I noticed (not sure if it does matter) is the difference in bytecode of the function solve:

Machine N.1 (NO ERROR)

> solve
function (a, b, ...) 
UseMethod("solve")
<bytecode: 0x208f570>
<environment: namespace:base>

Machine N.2 (WITH ERROR)

> solve
function (a, b, ...) 
UseMethod("solve")
<bytecode: 0x2c72718>
<environment: namespace:base>

The third point is: Machine N.1 has 2GB RAM and Machine N.2 (with error) has only 512MB RAM. Maybe it's a memory issue?

ets() optimizes the initial states, but traditional exponential smoothing does not.
HoltWinters() does not optimize the initial states, but it uses a non-standard heuristic initialization.
For teaching purposes, it would be good to have functions that did not optimize initial states but used standard heuristic initialization.

I envisage functions ses(), holt() and hw(), each with an argument "initial" taking values "optimal" or "simple". If initial=optimal, the function will just call ets(). Otherwise it will do something like HoltWinters() but with more standard initialization methods.

Mimic feature in arima.sim() for user to pass innovations.
Simple fix: look for call to rnorm and replace.

Hi Rob,

As per this stackoverflow thread http://stackoverflow.com/questions/20993650/dual-seasonal-cycles-in-ts-object/21007465?noredirect=1#comment31621171_21007465

I'm getting the following error when I try to use seasadj() with a bats object.

Error in seasadj(fit) : Object of unknown class

The error can be replicated (on my machine) with the following code:

x <-(21243,30908,39868,26176,23760,19531,16410,21444,14699,22149,18625,16410,14799,14397,20639,23156,17820,15605,14095,14095,12484)
x2 <- msts(x, seasonal.periods=c(3,7))
fit <- tbats(x2)
x.sa <- seasadj(fit)

See http://stats.stackexchange.com/a/55429/159

To compute forecasts and fitted values without re-estimating the model.

When attempting to run the main forecast function on a time series object exceptions are sometimes thrown depending on the size of the vector used to generate the time series object.

For example:

tsObj <- ts(tsVector[1:6], frequency=12, start=c(2013,4))

where tsVector is just a vector of values. Running this time series object through the forecast function like so:

pred <- forecast(tsObj, 1)

generates the following error:

Error in lsfit(1:maxn, y.sa[1:maxn]) : 0 responses, but only 2 variables
In addition: Warning message:
In lsfit(1:maxn, y.sa[1:maxn]) : 6 missing values deleted

Vectors of length 1,2,3 or >= 12 do not generate an error.

Also, just wanted to say that this package totally rocks! It has been instrumental in my work and I greatly appreciate it being made available to the wider world!

In ets you can't specify every trend model. The possibilities should include N, A, Ad, M, Md and Z and maybe Az and Mz (where the latter mean A or Ad and M or Md). If you specify the trend as A and damped = TRUE you can force Ad but if you specify damped = FALSE that does not force A. Also the need to specify two interacting arguments is confusing. Suggest they be combined and to retain that damped/undamped optimization capability perhaps Az and Mz as outlined.

I would like transfer function methodology to be added. We need to decide whether to use the Peter Young approach, as in the CAPTAIN toolbox in Matlab, or a state space approach as in Brockwell & Davis, or something else.

This code fails with the following output. Am I doing something wrong?
I guess it's an error in the underlying call to create the NN.

library(forecast)
library(tseries)
library(zoo)
dwa <- get.hist.quote(instrument="FAST", start="2009-12-31", end="2012-08-16")
returns <- dwa$Close / lag(dwa$Close, -1)
nn <- nnetar(coredata(returns))

Output:
Error in nnet.default(x, y, ...) : no weights to fit
Calls: nnetar -> avnnet -> -> nnet.default

See http://stackoverflow.com/questions/21722791/no-result-when-wrapping-working-forecast-statement-in-trycatch

--- rwf() fitted values contain NA where fit would be possible
data = 1:10
data[5]=NA
fitted(rwf(data))

[1] NA 1 2 3 NA NA 6 7 8 9

I would expect only one NA after the 3 in the output

--- auto.arima() 'simple' input gives an error
Some apparantly simple input cannot be modeled by arima. (Probably has to do with singular matrices? )

Fail scenario 1:
data=1:10
auto.arima(data)

Error in if (PVAL == min(tablep)) warning("p-value smaller than printed p-value") else warning("p-value greater than printed p-value") :
missing value where TRUE/FALSE needed

Fail scenario 2: (different error)
data=1:10
data[1]=data[1]*1.00000000000001
auto.arima(data)

Error in auto.arima(data) : No suitable ARIMA model found

Success scenario:
data=1:10
data[1]=data[1]*1.000000001
auto.arima(data)

--- auto.arima() with xreg gives error when data and xreg are too similar
Similar to previous, but now with regressors. Disturbing the target with a small pertubation fixes the problem.

Success scenario:
data = rep(1:6,6)
auto.arima(data)

Fail scenario:
data = rep(1:6,6)
datax = data
auto.arima(data,xreg=datax)

Success scenario:
data = rep(1:6,6)
datax = data
data[1]=data[1]*1.000000001
auto.arima(data,xreg=datax)

--- packageVersion("forecast") and loading the library give different numbers
packageVersion("forecast")
-> 4.6
library(forecast)
-> This is forecast 4.06

--- meanf gives different fit when input is timeseries and when input is vector
It seems that the one for the timeseries object is wrong:

fitted(meanf(rep(1:5,5)))

[1] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3

fitted(meanf(ts(rep(1:5,5),frequency=5)))

[1] NA 1.000000 1.500000 2.000000 2.500000 3.000000 2.666667 2.571429
[9] 2.625000 2.777778 3.000000 2.818182 2.750000 2.769231 2.857143 3.000000
[17] 2.875000 2.823529 2.833333 2.894737 3.000000 2.904762 2.863636 2.869565
[25] 2.916667

--- bats and tbats do not handle preceding missing values like ets
data=c(NA,1:10)
bats(data)

Error in if (any((y <= 0))) { : missing value where TRUE/FALSE needed

whereas with ets, successful with warning

data=c(NA,1:10)
ets(data)

Warning message:
In ets(data) :
Missing values encountered. Using longest contiguous portion of time series

Currently, these functions can handle xreg with method="arima", but ignore it when forecasting. Need to add a newxreg argument for forecasting purposes.