I am fairly new to statistics and R. I would like to know the process to determine the ARIMA parameters for my dataset. Can you help me figure out the same using R and theoretically (if possible)?
The data range from Jan-12 to Mar-14 and portrays the monthly sales. Here's the data set:
99 58 52 83 94 73 97 83 86 63 77 70 87 84 60 105 87 93 110 71 158 52 33 68 82 88 84
And here's the trend:
The data don't exhibit a trend, seasonal behavior or cyclicity.
In general, dig into an advanced time series analysis textbook (introductory books will usually direct you to just trust your software), like Time Series Analysis by Box, Jenkins & Reinsel. You may also find details on the Box-Jenkins procedure by googling. Note that there are other approaches than Box-Jenkins, e.g., AIC-based ones.
In R, you first convert your data into a ts (time series) object and tell R that the frequency is 12 (monthly data):
require(forecast)
sales <- ts(c(99, 58, 52, 83, 94, 73, 97, 83, 86, 63, 77, 70, 87, 84, 60, 105, 87, 93, 110, 71, 158, 52, 33, 68, 82, 88, 84),frequency=12)
You can plot the (partial) autocorrelation functions:
acf(sales)
pacf(sales)
These don't suggest any AR or MA behavior.
Then you fit a model and inspect it:
model <- auto.arima(sales)
model
See ?auto.arima for help. As we see, auto.arima chooses a simple (0,0,0) model, since it sees neither trend nor seasonality nor AR or MA in your data. Finally, you can forecast and plot the time series and forecast:
plot(forecast(model))
Look at ?forecast.Arima (note the capital A!).
This free online textbook is a great introduction to time series analysis and forecasting using R. Very much recommended.
Two things.Your time series is monthly,you need at least 4 years of data for a sensible ARIMA estimation, as reflected 27 points do not give the autocorrelation structure. This can also mean that your sales is affected by some external factors , rather than being correlated with its own value. Try to find out what factor affects your sales and is that factor being measured. Then you can run a regression or VAR (Vector Autoregression) to get forecasts.
If you absolutely don't have anything else other than these values , your best way is to use a exponential smoothing method to get a naive forecast. Exponential smoothing is available in R .
Secondly don't see the sales of a product in isolation, the sales of two products might be correlated for example increase in coffee sales can reflect decrease in tea sales. use the other product information to improve your forecast.
This typically happens with sales data in retail or supply chain. They dont show much of autocorrelation structure in the series. While on the other hand methods like ARIMA or GARCH typically work with stock market data or economic indices where you generally have autocorrelation.
This is really a comment but exceeds the allowable so I post it as a quasi-answer as it suggests the correct way to analyze time series data. .
The well-known fact but often ignored here and elsewhere is that the theoretical ACF/PACF which is used to formulate a tentative ARIMA model premises no Pulses/Level Shifts/Seasonal Pulses/Local Time Trends. Additionally it premises constant parameters and constant error variance over time. In this case the 21st observation (value=158) is easily flagged as an outlier/Pulse and a suggested adjustment of -80 yields a modified value of 78 . The resultant ACF/PACF of the modified series shows little or no evidence of stochastic (ARIMA) structure. In this case the operation was a success but the patient died. The sample ACF is based upon the covariance/variance and an unduly inflated/bloated variance yields a downwards bias to the ACF. Prof. Keith Ord once referred to this as the "Alice in Wonderland effect"
As it has been pointed out by Stephan Kolassa there is no much structure in your data. The autocorrelation functions do not suggest an ARMA structure (see acf(sales), pacf(sales)) and forecast::auto.arima does not choose any AR or MA order.
require(forecast)
require(tsoutliers)
fit1 <- auto.arima(sales, d=0, D=0, ic="bic")
fit1
#ARIMA(0,0,0) with non-zero mean
#Coefficients:
# intercept
# 81.3704
#s.e. 4.4070
Nevertheless, notice that the null of normality in the residuals is rejected at the 5% significance level.
JarqueBera.test(residuals(fit1))[[1]]
#X-squared = 12.9466, df = 2, p-value = 0.001544
Aside note: JarqueBera.test is based on function jarque.bera.test available in package tseries.
Including the additive outlier at observation 21 that is detected with tsoutliers renders normality in the residuals. Thus, the estimate of the intercept and the forecast are not affected by the outlying observation.
res <- tsoutliers::tso(sales, types=c("AO", "TC", "LS"),
args.tsmethod=list(ic="bic", d=0, D=0))
res
#ARIMA(0,0,0) with non-zero mean
#Coefficients:
# intercept AO21
# 78.4231 79.5769
#s.e. 3.3885 17.6072
#sigma^2 estimated as 298.5: log likelihood=-115.25
#AIC=236.49 AICc=237.54 BIC=240.38
#Outliers:
# type ind time coefhat tstat
#1 AO 21 2:09 79.58 4.52
JarqueBera.test(residuals(res$fit))[[1]]
#X-squared = 1.3555, df = 2, p-value = 0.5077
