How to fit a regression model with ARIMA errors on the seasonally adjusted component of a time series (in R)?

Question

I want to do these two things (combined) with a time series T:

forecast the seasonally adjusted component of T (STL used for the decomposition) and "add back" the seasonality (I assume that the seasonal component is unchanging, so I use naïve method for the seasonal component)
fit a regression model with ARIMA errors (exogenous regressors included in the formula)

In other words, I want to obtain forecasts using the seasonally adjusted component of T integrating an external predictor and "adding back" the seasonality.

I can do these two operations separately, but I can't get them to work in combination

Here is some toy examples:

First, load libraries and data:

library(forecast)
library(tsibble)
library(tibble)
library(tidyverse)
library(fable)
library(feasts)
library(fabletools)


us_change <- readr::read_csv("https://otexts.com/fpp3/extrafiles/us_change.csv") %>%
  mutate(Time = yearquarter(Time)) %>%
  as_tsibble(index = Time)

Example of fit and forecast with seasonally adjusted component of T:

model_def = decomposition_model(STL,
                                Consumption  ~ season(window = 'periodic') + trend(window = 13),
                                ARIMA(season_adjust ~ PDQ(0,0,0)),
                                SNAIVE(season_year),
                                dcmp_args = list(robust=TRUE)) 

fit <- us_change %>% model(model_def)

report(fit)

forecast(fit, h=8) %>% autoplot(us_change)

Example of regression model with ARIMA errors (Income as predictor):

model_def = ARIMA(Consumption ~ Income + PDQ(0,0,0))

fit <- us_change %>% model(model_def)

report(fit)

us_change_future <- new_data(us_change, 8) %>% mutate(Income = mean(us_change$Income))

forecast(fit, new_data = us_change_future) %>% autoplot(us_change)

These examples work, but I would like to do something like this:

model_def = decomposition_model(STL,
                                Consumption  ~ season(window = 'periodic') + trend(window = 13),
                                ARIMA(season_adjust ~ Income + PDQ(0,0,0)),
                                SNAIVE(season_year),
                                dcmp_args = list(robust=TRUE))


fit <- us_change %>% model(model_def)

report(fit)

us_change_future <- new_data(us_change, 8) %>% mutate(Income = mean(us_change$Income))

forecast(fit, new_data = us_change_future) %>% autoplot(us_change)

I get this output in the console:

> fit <- us_change %>% model(model_def)
Warning message:
1 error encountered for model_def
[1] object 'Income' not found

> 
> report(fit)
Series: Consumption 
Model: NULL model 
NULL model>

So I tried doing this in decomposition_model:

model_def = decomposition_model(STL,
                                Consumption  ~ season(window = 'periodic') + trend(window = 13),
                                ARIMA(season_adjust ~ us_change$Income + PDQ(0,0,0)),
                                SNAIVE(season_year),
                                dcmp_args = list(robust=TRUE))

No problem with the fit, but now I get an error in the forecast:

> forecast(fit, new_data = us_change_future) %>% autoplot(us_change)
Error in args_recycle(.l) : all(lengths == 1L | lengths == n) is not TRUE
In addition: Warning messages:
1: In cbind(xreg, intercept = intercept) :
  number of rows of result is not a multiple of vector length (arg 2)
2: In z[[1L]] + xm :
  longer object length is not a multiple of shorter object length

What am I doing wrong?

score 3 · Accepted Answer · answered Oct 08 '19 at 22:51

Nothing wrong with your code here, just something I hadn't considered people would do when making decomposition_model(). I've updated the decomposition modelling method to include exogenous regressors so that they can be used in component models (https://github.com/tidyverts/fabletools/commit/8dd505f6378327b8e93b8440ec17ecf9badf2561). If you update the package, your first attempt at modelling should work fine.

As for why the second attempt didn't work, the forecast method is finding us_change$Income and using that as the exogenous regressor for the future forecasts. This value has the length of us_change, which does not match the length of us_change_future, leading to the (confusing) error.

Reprex:

library(tidyverse)
library(tsibble)
library(fable)
library(feasts)

us_change <- readr::read_csv("https://otexts.com/fpp3/extrafiles/us_change.csv") %>%
  mutate(Time = yearquarter(Time)) %>%
  as_tsibble(index = Time)

model_def = decomposition_model(STL,
                                Consumption  ~ season(window = 'periodic') + trend(window = 13),
                                ARIMA(season_adjust ~ Income + PDQ(0,0,0)),
                                SNAIVE(season_year),
                                dcmp_args = list(robust=TRUE))


fit <- us_change %>% model(model_def)

report(fit)
#> Series: Consumption 
#> Model: STL decomposition model 
#> Combination: season_adjust + season_year
#> 
#> ========================================
#> 
#> Series: season_adjust 
#> Model: LM w/ ARIMA(1,0,2) errors 
#> 
#> Coefficients:
#>          ar1      ma1     ma2  Income  intercept
#>       0.6922  -0.5777  0.1975  0.2035     0.5993
#> s.e.  0.1163   0.1305  0.0755  0.0462     0.0883
#> 
#> sigma^2 estimated as 0.3234:  log likelihood=-157.39
#> AIC=326.77   AICc=327.24   BIC=346.16
#> 
#> Series: season_year 
#> Model: SNAIVE 
#> 
#> sigma^2: 0

us_change_future <- new_data(us_change, 8) %>% mutate(Income = mean(us_change$Income))

forecast(fit, new_data = us_change_future) %>% autoplot(us_change)

^{Created on 2019-10-09 by the reprex package (v0.2.1)}

How to fit a regression model with ARIMA errors on the seasonally adjusted component of a time series (in R)?

1 Answers1