Team CCAD Case Study 8: Crox stock ARIMA analysis
David Josephs, Andy Heroy, Carson Drake, Che’ Cobb
2020-03-02
library(quantmod)
library(tswge)
library(tswgewrapped)
library(dplyr)
library(magrittr)
library(tidyverse)
set.seed(666)
# this set of functions allows us to autoincrement our figures
counter <- function() {
x <- 0
return (
function() {
# Assigning outside of scope! The real purpose of <-
x <<- x+1
return(x)
}
)
}
# initialize a new counter
cnt <- counter()
# call the counter with a bolded figure caption
cap <- function(str) {
paste("**Figure", cnt(), ":**",str)
}1 Introduction
The purpose of this case study is to fit an ARIMA time series model to a stock of our choosing. Since we enjoy the finer consumer products available on the market, we chose to write our case study about CROX stock. Yes, that long forgotten footwear loved by all those who spend too much time on their feet in a given day.
Figure 1 : Crocs: utility and fashion for the business savvy
To accomplish this task, we’ve chosen to utilize the library (TSWGE) which was written and composed by Dr. Woodward, Gray, and Elliot, along the whamo brilliance of Dr. Bivin Sadler. This library contains many useful functions for time series analysis which we describe and explore in the analysis below.
2 Background
To begin, we’ll explain what an ARIMA model is comprised of, and how they are useful in time series forecasting techniques. An ARIMA model is comprised of two main components, one consisting of the autoregressive (AR) component, and one consisting of the moving average (MA) component. As we will see shortly, the integrated (I in ARIMA) as well as the seasonal portion are just special cases of the autoregressive model.
2.1 AR Models
The basic form of the AR part of ARIMA (with order \(1\)) is written as:
\[ X_t\sum_{j=0}^p\phi_jB^j = \epsilon_t \]
Where \(\mathrm{B}\) is the backshift operator such that \(\mathrm{B}^nX_t = X_{t-n}\), \(\phi_0 = 1\), and \(\epsilon_t \sim N(0, v)\). For our analysis, we can define the characteristic polynomial of an AR model in the following manner:
\[ \sum_{j=0}^p\phi_jZ^j = 0 \]
We can constrain this equation with the assumption that our AR processes are stationary. An AR process is stationary if and only if all the roots of its characteristic equation lie within the unit circle on the complex plane. That is to say for all roots of a characteristic AR polynomial, the AR process is stationary if for every root \(r\), \(\| \mathbf{r} \| < 1\). Visually this can be explained by viewing the unit circle on the complex plane:
Figure 2 : Unit circle in the complex plane, all complex roots must be WITHIN this circle
There is interesting behavior when the roots of the characteristic equation lie exactly on the unit circle. This is how autoregressive models capture integrated/wandering components and seasonal components.
2.1.1 The I in ARIMA
The Integrated (I) component of ARIMA occurs when the AR portion of the model has a linear root at exactly 1. This means that, using ARIMA(p,d,q) form, we can write a \(\left(0,1,0\right)\) model (\(d=1\)) as: \[ \left(1 - B\right)X_t = \epsilon_t \]
or more generally as
\[ \left(1-B\right)^dX_t = \epsilon_t \]
Where \(d\) is the integrated order of the model. It is important to note that in general, if you have \(d=2\), you are pushing it, while \(d \geq 3\) processes rarely occur on this earth. If you have \(d\) at a high value it is likely better to look to other models, for example high order AR models.
2.1.2 Seasonal ARIMA
Seasonal components also appear when the roots to the AR characteristic equations are at one, however instead of linear unit roots, Seasonal components have polynomial roots at one. That is, giving a seasonality \(s\):
\[ (1-B^s)X_t = \epsilon_t \]
2.2 MA models
While AR models are focused on the relationship between a realization \(X\) at time \(t\) and at time \(t-n\), MA models are focused on the noise portion of the equation. In general, MA models alone are not used for modeling real data, as they have some strange properties and are not quite as useful, but in conjunction with AR, integrated, and seasonal models you can build more powerful models of time series. An MA model is expressed just the same as an arima model, using polynomials:
\[ X_t = \epsilon_t\sum_{k=0}^q\theta_kB^k \]
Unlike the AR component, an MA model is a model of noise, which we will refer to as a General Linear Process (GLP). GLPs have a special property in that they are always stationary (in comparison to, for example, autoregressive processes, which are only stationary when their roots are bounded by the unit circle). However, MA processes have an undesirable characteristic: multiple processes (different sets of \(\theta\)) can have the same autocorrelation structure, that is they exhibit model multiplicity. That means we can have different models producing the same effects, which is not ideal. To get around this, we constrain the possible MA models, we must only use models with unique solutions. This characteristic is known as invertibility, and we can define the criteria for invertibility as having characteristic roots outside the unit circle.
2.3 General ARIMA models
To introduce the most general form of the ARIMA model, we will introduce some new notation for brevity:
\[\begin{align} \Phi_p(B) = \sum_{j=0}^p\phi_pB^p \\ \Theta_q(B) = \sum_{k=0}^q \theta_qB^q \end{align}\]This allows us to write seasonal ARIMA with order \((p,d,q)\) and seasonality \(s\) as:
\[ \Phi_p(B)(1-B)^d(1-B^s)X_t = \Theta_q(B)\epsilon_t \]
3 TSWGE
In this study we will be using the tswge library, which is paired with the book Applied Time Series Analysis with R[1], mentioned in the above paragraphs. The library contains a few main utilities. First, it has a series of functions for generating time series, gen.TYPEOFSERIES.wge. These allow for practice and developing intuition. Second, it has utilities for forecasting time series, fore.TYPEOFSERIES.wge. Finally, it has two irreplacable sets of utilities, those for estimating the order of a time series as well as exploring them(model identification), and those for estimating the coefficients once the model has been identified (parameter estimation). This library, while brilliant, is sometimes syntactically painful, so it is extended by the student made (and in active development!) tswgewrapped package (maintained by one of the authors of this report). We will use a mix of both libraries, and share source code for the functions used when appropriate, in order to explain exactly what we are doing.
4 Stock Analysis
4.1 Data Setup
As mentioned in our Introduction section, we are chiefly interested in the stock ticker CROX, for the company that makes crocs, everyone’s favorite footwear. Please note, due to the ongoing financial turmoil, we will be ignoring the last two weeks with coronavirus, as that is outside the scope of ARIMA. Lets first load in the data, and take a brief look at a candle chart. If you are actually considering investing in crocs, please note that they are going to be hit hard by the effects of the ongoing breakout, and read this news article
crox <- getSymbols("CROX",src = "yahoo")
cstock <- CROX
cstock_2years <-cstock['2018-02-01::2020-02-01']
candleChart(cstock_2years, theme = "white")
# assign to a data frame!
cstock_2years %<>% data.frame
Figure 3 : Candle stock of Croc’s stock over the last two years
Figure 4 : A meme for readers enjoyment
The first thing we must do is pick a forecast horizon. As we are forecasting stocks, it is unlikely and inappropriate to forecast stocks 15-20 days ahead. Therefore, we will focus on just one business week (5 days). To assess our model, we will use the following cross validation structure:
- cut the data at one forecast horizon (5 days), using the last 5 days as a holdout set
- Validate by using backtesting and getting the average square error(ASE). This technique is used extensively in Applied Time Series with R, and basically involves creating an ad hoc validation set. For more complex validation schemes, you could also use rolling windows, as discussed in Forecasting Principles and Practice (link in conclusion), but that is really unneccessary for the purposes of this study
- Assess models on the holdout set, using ASE
Lets go ahead and split our data, and isolate the closing price.
train_indices <- 1:(nrow(cstock_2years) - 5)
crox_train <- cstock_2years[train_indices,]$CROX.Close
crox_test <- cstock_2years[-train_indices,]$CROX.Close4.2 Exploratory Analysis
The first key part of doing any time series analysis is to look at the series, its ACF, and a periodigram (parzen window), which measures the frequency density of the series in decibels from 0 to the nyquist frequency (0.5):
x <- plotts.sample.wge(crox_train)
Figure 5 : Raw Data, ACF, on top, frequency periodigrams on bottom
What can we tell from this? We will work right to left. First, we do not see any clear seasonality in the line plot. This is expected. Crocs are a commodity which is in demand year round, so we do not expect any seasonalities. It is unclear as to whether or not there is an integrated or wandering trend, it could either have a low order integrated trend or a higher order AR term. The ACF tells the same story. These barely damping lines are highly indicative of integrated trends, but also of strong autoregressions (which if you recall, an integrated term can be viewed as an exceptionally strong autoregression). The parzen window on the bottom left tells the same story as the periodigram, but is in the authors’ opinion more digestible. This inverse curve shape is indicative of a wandering (integrated component), or of a higher order AR component (but more likely wandering trend).
We will continue our analysis as follows. We will consider two cases, the first that CROX stock is stationary, or that it stays the same over time. We will also consider the case that CROX is wandering (\(d>0\)). With these two setups, we will evaluate two techniques: visual identification of the model, and automated using an AIC and BIC based grid search technique, as favored by Woodward, Gray, and Elliot. First lets go ahead and difference the series, to see if that helps. We will use the artrans.wge function to difference it away
# define a function for a single order difference
diff1 <- . %>% artrans.wge(phi.tr = 1)
train_adj <- diff1(crox_train)
Figure 6 : First order ARIMA difference of CROX stock
par(mfrow=c(1,1))Here we see the series before and after the transformation. Lets talk about the ACF before the transformation first. This slowly damping ACF implies an AR model with real, positive roots, so we will look for that in our by hand analysis. The transformed ACF implies either an AR model or a simple ARMA model, or even simply noise, which would imply that CROX stock is simply performing a random walk. We will keep this in mind for our analysis.
4.3 Analysis of untransformed series
First, we will attempt to perform our analysis by hand. We said we wanted a high order AR model with positive, real roots. To confirm this, lets make a sample time series with a strong value of phi (a strongly autoregressive model). We will use the tswgewrapped generate function, which is a thin wrapper for the tswge version:
generate#> function (type, ...)
#> {
#> phrase <- paste0("tswge::gen.", rlang::enexpr(type), ".wge")
#> func <- rlang::parse_expr(phrase)
#> eval(rlang::call2(func, ...))
#> }
#> <bytecode: 0x000000001474e2a8>
#> <environment: namespace:tswgewrapped>acf(generate(arma, n=100, phi=(0.99), plot=FALSE), main="Untransformed ACF plot") 
Figure 7 : A strong autoregressive model looks fairly appropriate. In the wild, it would be a higher order AR model, but this is difficult to generate.
4.3.1 By hand analysis of the untransformed series
The first thing, looking at the ACF, which never goes negative, is we can definitely tell the AR process generating the data has strong, positive roots. This is apparent due to the lack of sign change and overall high value of ACF. Let us go ahead and try our first guess of model, and see if the roots match our assumptions. We will start with \(p = 3\). This is our choice because our generation of a random series with \(phi_1=0.99\) looked very similar to the autocorrelation of the original series. However, we would much prefer if the roots were a bit further away from the unit circle, while also exhibiting a stronger autoregressive effect. Therefore 3 is a good, simple start. We will first estimate the AR coefficients using the estimate function from tswgewrapped. The source code for the function is below:
estimate#> function (xs, p = NA, q = 0, type = "mle", ...)
#> {
#> if (all(is.na(p))) {
#> stop("You need to specify the AR order for the model.")
#> }
#> if (q > 0) {
#> return(tswge::est.arma.wge(xs, p, q, ...))
#> }
#> else {
#> return(tswge::est.ar.wge(xs, p, type, ...))
#> }
#> }
#> <bytecode: 0x0000000013a66e80>
#> <environment: namespace:tswgewrapped>Again, it is a thin wrapper for tswge, where it uses maximum likelihood estimation to quickly determine the ideal values for phi and theta.
train_hand_est <- estimate(crox_train, p=3)#>
#> Coefficients of Original polynomial:
#> 0.8452 0.1966 -0.0434
#>
#> Factor Roots Abs Recip System Freq
#> 1-0.9985B 1.0015 0.9985 0.0000
#> 1+0.2988B -3.3466 0.2988 0.5000
#> 1-0.1455B 6.8740 0.1455 0.0000
#>
#> This also conveniently prints out the roots of the polynomial. We see we have real roots, with two positive and one negative. The first root is rather close to the unit circle, which may not be desirable. Differencing is likely more appropriate, but we will continue for the sake of knowledge! We can also plot the ACF of this estimate, and see how it compares to our original data:
acf(generate(arma, n=length(crox_train), phi = train_hand_est$phi, plot=FALSE), main="ACF of sample plot")
Figure 8 : ACF of a sample series
This looks pretty good! We can check goodness of fit with two methods: first using residual plotting, and second using the ljung_box function, which performs the ljung box test. The ljung box test has a null hypothesis which the residuals are white noise, and an alternative that the residuals are correlated in some way. It is a portmanteau test. First we will look at the residuals visually. It is important to do this because the ljung box test, while the best of statistical tests for autocorrelation, has low statistical power.
plot_res <- function (res) {
par(mfrow = c(1, 2))
plot(res, type = "b")
title(main="Residual plot")
acf(res, main="ACF of Residuals")
par(mfrow = c(1, 1))
}
plot_res(train_hand_est$res)
Figure 9 : The residuals of the estimate. If it looks like noise, it is better. In this case, we still have those two extreme values.
This looks like an alright fit, however to double check we will go ahead with the portmanteau test. It is advisable to conduct this test at multiple lags, particularly 24 and 48 (suggested by box). Luckily, tswgewrapped has this all covered for us (source code included, no smoke and mirrors).
ljung_box#> function (x, p, q, k_val = c(24, 48))
#> {
#> ljung <- function(k) {
#> hush(tswge::ljung.wge(x = x, p = p, q = q, K = k))
#> }
#> sapply(k_val, ljung)
#> }
#> <bytecode: 0x000000001fb9c278>
#> <environment: namespace:tswgewrapped>ljung_box(train_hand_est$res, 3,0)#> [,1] [,2]
#> test "Ljung-Box test" "Ljung-Box test"
#> K 24 48
#> chi.square 34 61
#> df 21 45
#> pval 0.033 0.06So this is in agreement with our plot. Our fit is pretty close, however it is probably not a perfect fit. We did not model all the autocorrelation out of the data, thus this is not the best model. For a baseline, we will backtest and get the ASE of our model. For the sake of transparency, we will show the tswgewrapped source code below:
assess#> function (x, ...)
#> {
#> bcast <- fcst(x = x, ..., lastn = T)
#> ASE <- ase(x, bcast)
#> return(ASE)
#> }
#> <bytecode: 0x00000000200a3b88>
#> <environment: namespace:tswgewrapped>fcst#> function (type, ...)
#> {
#> phrase <- paste0("tswge::fore.", rlang::enexpr(type), ".wge")
#> func <- rlang::parse_expr(phrase)
#> eval(rlang::expr((!!func)(...)))
#> }
#> <bytecode: 0x0000000020156630>
#> <environment: namespace:tswgewrapped>ase#> function (x, xhat)
#> {
#> s <- length(x) - length(xhat$f) + 1
#> n <- length(x)
#> mean((xhat$f - x[s:n])^2)
#> }
#> <bytecode: 0x000000002022e4e0>
#> <environment: namespace:tswgewrapped>Again, it is just a thin wrapper for tswge, we forecast the series other than the last n.ahead points and get the ASE. Lets test it out:
ase_hand_un <- assess(crox_train, arma, n.ahead=5, phi = train_hand_est$phi)
title(main="5 step back forecast")
Figure 10 : Looks like a pretty good forecast
ase_hand_un#> [1] 1.1This is a pretty good forecast! Lets go ahead and automate one too:
4.3.2 Grid search analysis of the untransformed data
First, we will use AIC and BIC to determine the optimal values of p and q. Note the aic5.wge and tswgewrapped::aicbic are just criterion based grid searches. Lets go ahead and run it (note, we are using the tswgewrapped::hush function to avoid some annoying tswge output):
aics <- hush(aic5.wge(crox_train))
bics <- hush(aic5.wge(crox_train, type='bic'))
pander::pander(list(aics, bics))p q aic 11 3 1 -0.6896 10 3 0 -0.6871 14 4 1 -0.6858 13 4 0 -0.6853 16 5 0 -0.6849 p q bic 10 3 0 -0.6533 4 1 0 -0.6503 11 3 1 -0.6474 6 1 2 -0.6472 13 4 0 -0.643
It looks like we were pretty close! We will go ahead and try out an ARMA(3,1) model as well, for the sake of thoroughness. The analysis will be less verbose this time, as there is nothing new under the sun.
est_auto_un <- estimate(crox_train, p=3, q = 1)#>
#> Coefficients of Original polynomial:
#> 1.6557 -0.4727 -0.1834
#>
#> Factor Roots Abs Recip System Freq
#> 1-0.9979B 1.0021 0.9979 0.0000
#> 1-0.8692B 1.1504 0.8692 0.0000
#> 1+0.2114B -4.7306 0.2114 0.5000
#>
#> plot_res(est_auto_un$res)
Figure 11 : These residuals do not look much better, at a higher model complexity…
Lets go ahead and check with our Ljung Box too:
ljung_box(est_auto_un$res, 3, 1)#> [,1] [,2]
#> test "Ljung-Box test" "Ljung-Box test"
#> K 24 48
#> chi.square 31 57
#> df 20 44
#> pval 0.05 0.086This represents a very slight improvement, however at the risk of slightly overfitting. Lets see how much better our backtested forecast got:
ase_auto_un <- assess(crox_train, arma, n.ahead=5, phi = est_auto_un$phi, theta = est_auto_un$theta)
title(main="5 step back forecast")
Figure 12 : Auto generated backtesting
ase_auto_un#> [1] 1.3Interestingly, this model did slightly worse in backtesting. It is likely that the MA term is not totally appropriate with the series, as we did not see any indicators in the plots or anywhere else in the data that suggested MA, and we are just modeling noise. We will keep it for comparison. Lets move on to the differenced series now.
4.4 Analysis of the differenced series
We will now look at the first order differenced data in more detail. First lets look at the sample plot:
# assignment to save the grader from tswge output
x <- plotts.sample.wge(train_adj)
Figure 13 : The differenced series is suggestive of MA or AR with complex conjugate and/or negative roots
As the caption states, the oscillatory nature of the ACF, suggests AR with either conjugate or negative roots. What is more interesting is the parzen window on the bottom left. This sharply lumpy state can be indicative of a couple of things: first, it is possible we are left with just white noise and not enough data to flatten out the parzen window, and second: a mixture of AR and MA components. AR with complex conjugate roots would make the peaks, while MA would make the troughs. Since conjugate roots come in pairs, and the oscillations continue for a long time, we can assume both the AR and the MA components have an even number. The autocorrelations continue to extend and oscillate for a long period of time, with a slow sinusoidal damping pattern. This suggests to us that the AR component outweighs the MA component (and in general we prefer more AR than MA, as pure AR models are in general simpler and should generalize better). We will go with ARMA(4,2) because of these reasons (the parzen window is such a clear sign of an ARMA model).
est_hand_diff <- estimate(train_adj, p=4,q=2)#>
#> Coefficients of Original polynomial:
#> -0.8924 -0.9981 -0.0707 0.0710
#>
#> Factor Roots Abs Recip System Freq
#> 1+0.7623B+0.9720B^2 -0.3921+-0.9354i 0.9859 0.3132
#> 1+0.3431B -2.9147 0.3431 0.5000
#> 1-0.2131B 4.6937 0.2131 0.0000
#>
#> These roots make a lot of sense. First of all, we have the complex conjugate, which explains the peak in the parzen window at 0.3132. The second root is negative, with a peak near the nyquist frequency in the parzen window. The third root is relatively weak, and shows the remaining wandering component to the series. We can also check the roots of the MA component of the model using the brilliant factor.wge function (note this doesn’t explain the interaction between the two sets of roots, but it definitely helps)
factor.wge(est_hand_diff$theta)#>
#> Coefficients of Original polynomial:
#> -0.7350 -0.9402
#>
#> Factor Roots Abs Recip System Freq
#> 1+0.7350B+0.9402B^2 -0.3909+-0.9544i 0.9697 0.3119
#>
#> We see another complex conjugate root, with a frequency of 0.3119. This is slightly less than the peak from the AR part, which probably explains the sharpness of the central peak. Again, we cannot model the interplay between the two sets of roots with simple factoring. Lets go ahead and run our tests:
plot_res(est_hand_diff$res)
Figure 14 : Not much change…
ljung_box(est_hand_diff$res, p=4,q=2)#> [,1] [,2]
#> test "Ljung-Box test" "Ljung-Box test"
#> K 24 48
#> chi.square 25 47
#> df 18 42
#> pval 0.13 0.29Our residual plot looks about the same, however finally we can safely accept the null hypothesis and say our residuals are stationary! Lets test the forecast out:
ase_hand_diff <- assess(crox_train, aruma, n.ahead=5, phi = est_hand_diff$phi, theta =est_hand_diff$theta, d=1)
title(main="Backtesting ARIMA")
Figure 15 : Worst forecast so far!
ase_hand_diff#> [1] 1.5This is not a brilliant forecast, relative to what we have done previously. This is largely in part due to some of the facts regarding ARIMA. The integrated component of ARIMA is more or less a random walk, and the expected value of a random walk is just wherever you started. So, a forecast with d=1 and nothing else is just going to predict the last value over and over and over. When we add in the ARMA components, we are just adding a little wiggle to that same pattern. This means even if it is the most appropriate model, ARIMA is not always the best forecast. We will investigate this further in the conclusion section. Lets now try to automate the choice of model using information criterion again.
4.4.1 Automated Analysis of differenced series:
aics <- hush(aic5.wge(train_adj))
bics <- hush(aic5.wge(train_adj, type='bic'))
pander::pander(list(aics, bics))p q aic 15 4 2 -0.7054 4 1 0 -0.7002 3 0 2 -0.6984 2 0 1 -0.6974 7 2 0 -0.6972 p q bic 4 1 0 -0.6832 2 0 1 -0.6804 3 0 2 -0.673 7 2 0 -0.6718 5 1 1 -0.6715
It looks like a simple ARMA(1,0) model is most appropriate (other than our 4,2 model, which did not fare so well). I like the idea of this simple model, lets estimate the parameters and check out the residuals
est_auto_diff <- estimate(train_adj, p =1)#>
#> Coefficients of Original polynomial:
#> -0.1690
#>
#> Factor Roots Abs Recip System Freq
#> 1+0.1690B -5.9183 0.1690 0.5000
#>
#> plot_res(est_auto_diff$res)
Figure 16 : Final residual plot looks slightly worse than the previous
ljung_box(est_auto_diff$res, p=1, q=0)#> [,1] [,2]
#> test "Ljung-Box test" "Ljung-Box test"
#> K 24 48
#> chi.square 34 58
#> df 23 47
#> pval 0.072 0.12This got us slightly worse than the previous model as far as modeling the noise, but these portmanteau tests are simply guidelines, as they are pretty weak. We can proceed with assessing the model:
ase_auto_diff <- assess(crox_train, aruma, n.ahead=5, phi = est_auto_diff$phi, d=1)
title(main="Backtesting ARIMA")
Figure 17 : Final backtest, looks like a classic d=1 model
ase_auto_diff#> [1] 1.3This model performed all right. It is simpler than the previous model, and performed better than it, so I am more likely to choose it over the 4,2 model. Lets go ahead and test all the models on the holdout set:
5 Conclusions
We will start by testing the first model made, by hand using no seasonal differencing on the test set:
f_un_h <- fcst( arma, crox_train, n.ahead=5, phi = train_hand_est$phi)
lines(x=(length(train_indices)):nrow(cstock_2years), y=c(crox_train[length(crox_train)],crox_test), col='red')
title(main="Hand made forecast")
Figure 18 : An alright forecast
ase(crox_test, f_un_h)#> [1] 3.3This is not a brilliant forecast, but it is not horrible, we have to keep in mind that at this point in time there is already worry over a global pandemic. Lets look at the automatically made forecast and then discuss ARMA forecasts in general:
f_a_h <- fcst(arma,crox_train, n.ahead=5, phi = est_auto_un$phi, theta = est_auto_un$theta)
lines(x=(length(train_indices)):nrow(cstock_2years), y=c(crox_train[length(crox_train)],crox_test), col='red')
title(main="Grid/Criterion made forecast")
Figure 19 : A better forecast!
ase(crox_test, f_a_h)#> [1] 3.1This did slightly better! We were able to capture more of the downward trend more quickly with the stronger fit by including the MA term. Lets briefly discuss the other advantages of not including the I term. With ARMA models, in general we are suggesting that things are going to trend towards the mean, as they are stationary models. For goods which do not explode in demand, but are still pretty constantly demanded, such as crocs, this is a reasonable assumption. We can display this by making a very long forecast with our ARMA model:
f_l <- fcst(arma,crox_train, n.ahead=1600, phi = est_auto_un$phi, theta = est_auto_un$theta)
abline(h=mean(crox_train), col='red')
title(main="Extreme ARMA forecasting")
Figure 20 : Trending towards the mean(in red)…
As t approaches infinity, we will eventually just start to forecast the mean over and over. So ARMA models are conservative, saying things will stay the same as they have in the past. Lets check out our ARIMA models now.
f_d_h <- fcst(aruma, crox_train, n.ahead=5, phi = est_hand_diff$phi, theta =est_hand_diff$theta, d=1)
lines(x=(length(train_indices)):nrow(cstock_2years), y=c(crox_train[length(crox_train)],crox_test), col='red')
title(main="Hand made first order difference forecast")
Figure 21 : The slope appears way off
ase(crox_test, f_d_h)#> [1] 3.5Our worst forecast so far. We will also go ahead and check out the automated one:
f_a_d <- fcst(aruma, crox_train, n.ahead=5, phi = est_auto_diff$phi, d=1)
lines(x=(length(train_indices)):nrow(cstock_2years), y=c(crox_train[length(crox_train)],crox_test), col='red')
title(main="BIC made first order difference forecast")
Figure 22 : Our dead simple model had an even worse forecast…
ase(crox_test, f_a_d)#> [1] 3.8This is the worst model we made. Why did the differenced models perform so poorly? Isn’t ARIMA a powerful tool?? The reason is this: When you make an ARIMA forecast, you are trying to predict a random walk, or an AR process at with a root at 1. When you do this, the math happens to fall in place such that the model just predicts the last value seen repeatedly (with d=2, the slope of the last 2 values, with d=3, you are doing something wrong). We can test this by making another extremely long forecast:
f_l <- fcst(aruma, crox_train, n.ahead=1600, phi = est_auto_diff$phi, d=1)
abline(h=mean(crox_train), col='red')
title(main="Extreme ARIMA forecasting")
Figure 23 : I predict things will be exactly the same as they are right now for the rest of time
ARIMA is a powerful tool, and incorporating the random walk/integrated component can help us build powerful models (especially with seasonality involved), however when the only thing that seems to be going on is an integrated component (or strong AR), it is often better to underfit that as a strong AR model than overfit it as an ARIMA model. All in all, the toolset given by tswge is a powerful one. For further reference, please refer to the study guide written by the authors of this report (link), the tswgewrapped documnetation on rdrr.io, the Applied Time Series Analysis with R book, and the brilliant book by Rob J Hyndman, the bible for forecasting, fpp2.
5.1 References
[1] Wayne A. Woodward, Henry L. Gray, Alan C. Elliott, Applied Time Series Analysis in R. Boca Raton: CRC Press.