Solved – forecasting multivariate time series (with categorical variables) in R

categorical dataforecastingmultivariate analysisrtime series

I want to forecast future(next 20 days) sales with sample dataset.

This is just a sample data and the actual data is from Jan 2014 to Dec 2016. As you can see, sales tend to increase as time goes by, and usually it gets higher when variable(c) is 'no'.

Time series forecasting cases that were posted on Stackoverflow or other Websites were either univariate time series or multivariate time series without categorical variables. I definitely don't want to ignore the variable(c) in this case.

I tried to use Random forest and XGBOOST but the result was really bad. It is because those decision-tree based models don't work well for extrapolation analysis. I am pretty sure regression based model will work well for this case but linear regression can't read the categorical variable and also, it won't capture the seasonality and weekdays trend.

To sum up, it would be great if you can recommend me a good future forecasting model that satisfies those 3 things:

Captures the sales increasing trend
Captures the seasonality and weekdays trend
Treats not only numerical variables but also categorical variables

Thank you.

Best Answer

A linear regression function like lm in R will correctly interpret your categorical variable. Since it ony has two states, you might recode it as 0/1 anyway.

I think what you want is an ARIMAX model (ARIMA with eXogeneous variables). Alternatively, you might want to use a state-space model which also allows to introduce regressors (even with time-varying coefficients). For the first approach, if you are using R you might turn to package forecast (Arima function). For the second, you might look at packages dlm or KFAS.

Related Solutions

Solved – stochastic vs deterministic trend/seasonality in time series forecasting

1) As regards your first question, some tests statistics have been developed and discussed in the literature to test the null of stationarity and the null of a unit root. Some of the many papers that were written on this issue are the following:

Related to the trend:

Dickey, D. y Fuller, W. (1979a), Distribution of the estimators for autoregressive time series with a unit root, Journal of the American Statistical Association 74, 427-31.
Dickey, D. y Fuller, W. (1981), Likelihood ratio statistics for autoregressive time series with a unit root, Econometrica 49, 1057-1071.
Kwiatkowski, D., Phillips, P., Schmidt, P. y Shin, Y. (1992), Testing the null hypothesis of stationarity against the alternative of a unit root: How sure are we that economic time series have a unit root?, Journal of Econometrics 54, 159-178.
Phillips, P. y Perron, P. (1988), Testing for a unit root in time series regression, Biometrika 75, 335-46.
Durlauf, S. y Phillips, P. (1988), Trends versus random walks in time series analysis, Econometrica 56, 1333-54.

Related to the seasonal component:

Hylleberg, S., Engle, R., Granger, C. y Yoo, B. (1990), Seasonal integration and cointegration, Journal of Econometrics 44, 215-38.
Canova, F. y Hansen, B. E. (1995), Are seasonal patterns constant over time? a test for seasonal stability, Journal of Business and Economic Statistics 13, 237-252.
Franses, P. (1990), Testing for seasonal unit roots in monthly data, Technical Report 9032, Econometric Institute.
Ghysels, E., Lee, H. y Noh, J. (1994), Testing for unit roots in seasonal time series. some theoretical extensions and a monte carlo investigation, Journal of Econometrics 62, 415-442.

The textbook Banerjee, A., Dolado, J., Galbraith, J. y Hendry, D. (1993), Co-Integration,Error Correction, and the econometric analysis of non-stationary data, Advanced Texts in Econometrics. Oxford University Press is also a good reference.

2) Your second concern is justified by the literature. If there is a unit root test then the traditional t-statistic that you would apply on a linear trend does not follow the standard distribution. See for example, Phillips, P. (1987), Time series regression with unit root, Econometrica 55(2), 277-301.

If a unit root exists and is ignored, then the probability of rejecting the null that the coefficient of a linear trend is zero is reduced. That is, we would end up modelling a deterministic linear trend too often for a given significance level. In the presence of a unit root we should instead transform the data by taking regular differences to the data.

3) For illustration, if you use R you can do the following analysis with your data.

x <- structure(c(7657, 5451, 10883, 9554, 9519, 10047, 10663, 10864, 
  11447, 12710, 15169, 16205, 14507, 15400, 16800, 19000, 20198, 
  18573, 19375, 21032, 23250, 25219, 28549, 29759, 28262, 28506, 
  33885, 34776, 35347, 34628, 33043, 30214, 31013, 31496, 34115, 
  33433, 34198, 35863, 37789, 34561, 36434, 34371, 33307, 33295, 
  36514, 36593, 38311, 42773, 45000, 46000, 42000, 47000, 47500, 
  48000, 48500, 47000, 48900), .Tsp = c(1, 57, 1), class = "ts")

First, you can apply the Dickey-Fuller test for the null of a unit root:

require(tseries)
adf.test(x, alternative = "explosive")
#   Augmented Dickey-Fuller Test
#   Dickey-Fuller = -2.0685, Lag order = 3, p-value = 0.453
#   alternative hypothesis: explosive

and the KPSS test for the reverse null hypothesis, stationarity against the alternative of stationarity around a linear trend:

kpss.test(x, null = "Trend", lshort = TRUE)
#   KPSS Test for Trend Stationarity
#   KPSS Trend = 0.2691, Truncation lag parameter = 1, p-value = 0.01

Results: ADF test, at the 5% significance level a unit root is not rejected; KPSS test, the null of stationarity is rejected in favour of a model with a linear trend.

Aside note: using lshort=FALSE the null of the KPSS test is not rejected at the 5% level, however, it selects 5 lags; a further inspection not shown here suggested that choosing 1-3 lags is appropriate for the data and leads to reject the null hypothesis.

In principle, we should guide ourselves by the test for which we were able to the reject the null hypothesis (rather than by the test for which we did not reject (we accepted) the null). However, a regression of the original series on a linear trend turns out to be not reliable. On the one hand, the R-square is high (over 90%) which is pointed in the literature as an indicator of spurious regression.

fit <- lm(x ~ 1 + poly(c(time(x))))
summary(fit)
#Coefficients:
#                 Estimate Std. Error t value Pr(>|t|)    
#(Intercept)       28499.3      381.6   74.69   <2e-16 ***
#poly(c(time(x)))  91387.5     2880.9   31.72   <2e-16 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#Residual standard error: 2881 on 55 degrees of freedom
#Multiple R-squared:  0.9482,   Adjusted R-squared:  0.9472 
#F-statistic:  1006 on 1 and 55 DF,  p-value: < 2.2e-16

On the other hand, the residuals are autocorrelated:

acf(residuals(fit)) # not displayed to save space

Moreover, the null of a unit root in the residuals cannot be rejected.

adf.test(residuals(fit))
#   Augmented Dickey-Fuller Test
#Dickey-Fuller = -2.0685, Lag order = 3, p-value = 0.547
#alternative hypothesis: stationary

At this point, you can choose a model to be used to obtain forecasts. For example, forecasts based on a structural time series model and on an ARIMA model can be obtained as follows.

# StructTS
fit1 <- StructTS(x, type = "trend")
fit1
#Variances:
# level    slope  epsilon  
#2982955        0   487180 
# 
# forecasts
p1 <- predict(fit1, 10, main = "Local trend model")
p1$pred
# [1] 49466.53 50150.56 50834.59 51518.62 52202.65 52886.68 53570.70 54254.73
# [9] 54938.76 55622.79

# ARIMA
require(forecast)
fit2 <- auto.arima(x, ic="bic", allowdrift = TRUE)
fit2
#ARIMA(0,1,0) with drift         
#Coefficients:
#         drift
#      736.4821
#s.e.  267.0055
#sigma^2 estimated as 3992341:  log likelihood=-495.54
#AIC=995.09   AICc=995.31   BIC=999.14
#
# forecasts
p2 <- forecast(fit2, 10, main = "ARIMA model")
p2$mean
# [1] 49636.48 50372.96 51109.45 51845.93 52582.41 53318.89 54055.37 54791.86
# [9] 55528.34 56264.82

A plot of the forecasts:

par(mfrow = c(2, 1), mar = c(2.5,2.2,2,2))
plot((cbind(x, p1$pred)), plot.type = "single", type = "n", 
  ylim = range(c(x, p1$pred + 1.96 * p1$se)), main = "Local trend model")
grid()
lines(x)
lines(p1$pred, col = "blue")
lines(p1$pred + 1.96 * p1$se, col = "red", lty = 2)
lines(p1$pred - 1.96 * p1$se, col = "red", lty = 2)
legend("topleft", legend = c("forecasts", "95% confidence interval"), 
  lty = c(1,2), col = c("blue", "red"), bty = "n")
plot((cbind(x, p2$mean)), plot.type = "single", type = "n", 
  ylim = range(c(x, p2$upper)), main = "ARIMA (0,1,0) with drift")
grid()
lines(x)
lines(p2$mean, col = "blue")
lines(ts(p2$lower[,2], start = end(x)[1] + 1), col = "red", lty = 2)
lines(ts(p2$upper[,2], start = end(x)[1] + 1), col = "red", lty = 2)

trend forecasts

The forecasts are similar in both cases and look reasonable. Notice that the forecasts follow a relatively deterministic pattern similar to a linear trend, but we did not modelled explicitly a linear trend. The reason is the following: i) in the local trend model, the variance of the slope component is estimated as zero. This turns the trend component into a drift that has the effect of a linear trend. ii) ARIMA(0,1,1), a model with a drift is selected in a model for the differenced series.The effect of the constant term on a differenced series is a linear trend. This is discussed in this post.

You may check that if a local model or an ARIMA(0,1,0) without drift are chosen, then the forecasts are a straight horizontal line and, hence, would have no resemblance with the observed dynamic of the data. Well, this is part of the puzzle of unit root tests and deterministic components.

Edit 1 (inspection of residuals): The autocorrelation and partial ACF do not suggest a structure in the residuals.

resid1 <- residuals(fit1)
resid2 <- residuals(fit2)
par(mfrow = c(2, 2))
acf(resid1, lag.max = 20, main = "ACF residuals. Local trend model")
pacf(resid1, lag.max = 20, main = "PACF residuals. Local trend model")
acf(resid2, lag.max = 20, main = "ACF residuals. ARIMA(0,1,0) with drift")
pacf(resid2, lag.max = 20, main = "PACF residuals. ARIMA(0,1,0) with drift")

ACF-PACF

As IrishStat suggested, checking for the presence of outliers is also advisable. Two additive outliers are detected using the package tsoutliers.

require(tsoutliers)
resol <- tsoutliers(x, types = c("AO", "LS", "TC"), 
  remove.method = "bottom-up", 
  args.tsmethod = list(ic="bic", allowdrift=TRUE))
resol
#ARIMA(0,1,0) with drift         
#Coefficients:
#         drift        AO2       AO51
#      736.4821  -3819.000  -4500.000
#s.e.  220.6171   1167.396   1167.397
#sigma^2 estimated as 2725622:  log likelihood=-485.05
#AIC=978.09   AICc=978.88   BIC=986.2
#Outliers:
#  type ind time coefhat  tstat
#1   AO   2    2   -3819 -3.271
#2   AO  51   51   -4500 -3.855

Looking at the ACF, we can say that, at the 5% significance level, the residuals are random in this model as well.

par(mfrow = c(2, 1))
acf(residuals(resol$fit), lag.max = 20, main = "ACF residuals. ARIMA with additive outliers")
pacf(residuals(resol$fit), lag.max = 20, main = "PACF residuals. ARIMA with additive outliers")

enter image description here

In this case, the presence of potential outliers does not appear to distort the performance of the models. This is supported by the Jarque-Bera test for normality; the null of normality in the residuals from the initial models (fit1, fit2) is not rejected at the 5% significance level.

jarque.bera.test(resid1)[[1]]
# X-squared = 0.3221, df = 2, p-value = 0.8513
jarque.bera.test(resid2)[[1]]
#X-squared = 0.426, df = 2, p-value = 0.8082

Edit 2 (plot of residuals and their values) This is how the residuals look like:

residuals

And these are their values in a csv format:

0;6.9205
-0.9571;-2942.4821
2.6108;4695.5179
-0.5453;-2065.4821
-0.2026;-771.4821
0.1242;-208.4821
0.1909;-120.4821
-0.0179;-535.4821
0.1449;-153.4821
0.484;526.5179
1.0748;1722.5179
0.3818;299.5179
-1.061;-2434.4821
0.0996;156.5179
0.4805;663.5179
0.8969;1463.5179
0.4111;461.5179
-1.0595;-2361.4821
0.0098;65.5179
0.5605;920.5179
0.8835;1481.5179
0.7669;1232.5179
1.4024;2593.5179
0.3785;473.5179
-1.1032;-2233.4821
-0.3813;-492.4821
2.2745;4642.5179
0.2935;154.5179
-0.1138;-165.4821
-0.8035;-1455.4821
-1.2982;-2321.4821
-1.9463;-3565.4821
-0.1648;62.5179
-0.1022;-253.4821
0.9755;1882.5179
-0.5662;-1418.4821
-0.0176;28.5179
0.5;928.5179
0.6831;1189.5179
-1.8889;-3964.4821
0.3896;1136.5179
-1.3113;-2799.4821
-0.9934;-1800.4821
-0.4085;-748.4821
1.2902;2482.5179
-0.0996;-657.4821
0.5539;981.5179
2.0007;3725.5179
1.0227;1490.5179
0.27;263.5179
-2.336;-4736.4821
1.8994;4263.5179
0.1301;-236.4821
-0.0892;-236.4821
-0.1148;-236.4821
-1.1207;-2236.4821
0.4801;1163.5179

Solved – Multilinear regression vs. Time Series

The most popular methods for time series forecasting are ARIMA and Holt-Winters. Holt-Winters is simpler and computationally less expensive. It is especially designed for seasonality. If you use ARIMA you must use seasonal ARIMA. Some studies show that, in real life data, Holt-Winters and seasonal ARIMA not only have a similar accuracy but predict similar values (http://webarchive.nationalarchives.gov.uk/20160105160709/http://www.ons.gov.uk/ons/guide-method/ukcemga/ukcemga-publications/publications/archive/from-holt-winters-to-arima-modelling--measuring-the-impact-on-forecasting-errors-for-components-of-quarterly-estimates-of-public-service-output.pdf).

Your problem is multi-variate time series, not just $y(t)$ but $y(X,t)$.

An idea is to try something like a linear regression with time varying coefficients. Think of a linear model without time : $y(X)=\beta X$ where $\beta$ is a vector. Then introduce time : $y(X,t)=\beta(t) X$. Then see $\beta(t)$ as a vector time series.

But you can't observe $\beta(t)$.

I have not read any article addressing this problem even though it is somehow a problem of growing importance : studying time and other parameters together.

If at each time you have a lot of (y,X) samples covering well the possible values of X, then you could possibly estimate $\beta(t)$ at every time and then apply a time series method to this vector. But I guess you data is not like this.

Added:

The ultimate solution where you can really work on $\beta(t)$ as a hidden parameter is the Kalman filter. If $\beta(t)$ is just a random walk, then the state of Kalman filter is just $\beta(t)$ and this would be not so difficult to implement, using observations to update the filter. Next step would be to use an ARIMA or HW model inside the filter but this starts being hardcore maths. I don't know of anything like this packaged as an easy-to-use tool. You can read more about it: What are disadvantages of state-space models and Kalman Filter for time-series modelling?

Best Answer

Related Solutions

Solved – stochastic vs deterministic trend/seasonality in time series forecasting

Solved – Multilinear regression vs. Time Series

Related Question