Time Series and ARMA models

A time series is a series of data points indexed (or listed or graphed) in time order. Most commonly, a time series is a sequence taken at successive equally spaced points in time. Thus it is a sequence of discrete-time data. Examples of time series can be taken from our real life such as are weather forecast, monthly grocery sales, website traffic, stock market trend and predictions and many more.

Time series analysis comprises methods for analyzing time series data in order to extract meaningful statistics and other characteristics of the data. Time series forecasting is the use of a model to predict future values based on previously observed values. In lay-man’s language forecasting is how time series like this will continue same way in future. Whilst there are many models that can be applied to time series, for the purpose of this vignette we will focus on most commonly used models - Auto-Regression(AR) and Moving Average(MA) model.

Stationarity

Before we start applying the beautiful pre-made functions made available to us in R, we should understand the most important charcteristic of time series. Why do we care about ‘stationarity’ of a time series?

The reason this is brought to attention is that: until unless the time series is stationary, we cannot build a time series model! specifically, AR or MA are not applicable on non-stationary series.There are three basic criterion for a series to be classified as stationary series :

  1. The mean of the series should not be a function of time rather should be a constant.
  2. The variance of the series should not a be a function of time. This property is known as homoscedasticity.
  3. The covariance of the i th term and the (i + m) th term should not be a function of time.

In cases where the stationary criterion are violated, the first requisite becomes to stationarize the time series and then try stochastic models to predict this time series. There are multiple ways of bringing this stationarity - three commonly used are: 1. Detrending : Here, we simply remove the trend component from the time series by using log() 2. Differencing : This is the commonly used technique to remove non-stationarity. Here we try to model the differences of the terms and not the actual term. For this we use funtion diff() 3. Seasonality : Seasonality can easily be incorporated in the ARIMA model directly, with argument ’S’.

Introduction to ARMA Time Series Modeling

AR - Auto Regressive:

The autoregressive (AR) model is arguably the most widely used time series model. It shares the very familiar interpretation of a simple linear regression, but here each observation is regressed on the previous observation.

In simple language - todays observation is regressed on yesterdays observation. 

equation for AR:  x(t) = alpha * x(t – 1) + error (t)

This equation is known as AR(1) formulation. The numeral one (1) denotes that the next instance is solely dependent on the previous instance. The alpha is a coefficient which we seek so as to minimize the error function. Notice that x(t- 1) is indeed linked to x(t-2) in the same fashion. Hence, any shock to x(t) will gradually fade off in future.

  1. If alpha=0 then x(t) is simply white noise
  2. Else x(t) is auto correlated and not stationary
  3. Large alpha lead to greater autocorrelation
  4. Negative values of alpha result in oscillatory time series

MA - Moving Averages:

The simple moving average (MA) model is a parsimonious time series model used to account for very short-run autocorrelation. It does have a regression like form, but here each observation is regressed on the previous innovation, which is not actually observed. A weighted sum of previous and current noise is called Moving Average (MA) model.

In simple language - SMA is todays observation is regressed on yesterdays noise.

equation for MA: x(t) = beta * error(t-1) + error (t)

This equation is known as MA(1) formulation. The numeral one (1) denotes that the next instance is solely dependent on the previous instance. The beta is a coefficient which we seek so as to minimize the error function.

  1. If beta is 0 then x(t) is simply a white noise process
  2. Else x(t) is auto correlated
  3. Large values of beta lead to larger correlation and
  4. Negative values of beta result in Oscillatory time series
  5. Auto correlation on MA is only for 1 lag

Difference between AR and MA models: The primary difference between an AR and MA model is based on the correlation between time series objects at different time points. The correlation between x(t) and x(t-n) for n > order of MA is always zero. This directly flows from the fact that covariance between x(t) and x(t-n) is zero for MA models (something which we refer from the example taken in the previous section). However, the correlation of x(t) and x(t-n) gradually declines with n becoming larger in the AR model. This difference gets exploited irrespective of having the AR model or MA model. The correlation plot can give us the order of MA model.

Exploration of Time Series Data in R

Exploring data becomes most important in a time series model – without this exploration, you will not know whether a series is stationary or not.

To elaborate the ARMA models we will use an inbuilt data set of R called AirPassengers. The dataset consists of monthly totals of international airline passengers, 1949 to 1960.

STEP 1:

Load the data. Then the first pre-requisite is, no prizes for guessing - the data should be a time series data and to check that you can use function ‘is.ts()’.

#Loading the Data Set
data("AirPassengers")

#This tells you that the data series is in a time series format
is.ts(AirPassengers)
## [1] TRUE

STEP 2:

Now that we know that the data is time series we should do some data exploration. Functions print() and summary() are used to get the overview of the data. The start() and end() functions return the time index of the first and last observations, respectively. The time() function calculates a vector of time indices, with one element for each time index on which the series was observed. Finally, the frequency() function returns the number of observations per unit time.

#This will give us the structure of our data
print(AirPassengers)
##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 112 118 132 129 121 135 148 148 136 119 104 118
## 1950 115 126 141 135 125 149 170 170 158 133 114 140
## 1951 145 150 178 163 172 178 199 199 184 162 146 166
## 1952 171 180 193 181 183 218 230 242 209 191 172 194
## 1953 196 196 236 235 229 243 264 272 237 211 180 201
## 1954 204 188 235 227 234 264 302 293 259 229 203 229
## 1955 242 233 267 269 270 315 364 347 312 274 237 278
## 1956 284 277 317 313 318 374 413 405 355 306 271 306
## 1957 315 301 356 348 355 422 465 467 404 347 305 336
## 1958 340 318 362 348 363 435 491 505 404 359 310 337
## 1959 360 342 406 396 420 472 548 559 463 407 362 405
## 1960 417 391 419 461 472 535 622 606 508 461 390 432
#This will give us summary of our data
summary(AirPassengers)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   104.0   180.0   265.5   280.3   360.5   622.0
#Starting index, end index
start(AirPassengers)
## [1] 1949    1
end(AirPassengers)
## [1] 1960   12
time(AirPassengers)
##           Jan      Feb      Mar      Apr      May      Jun      Jul
## 1949 1949.000 1949.083 1949.167 1949.250 1949.333 1949.417 1949.500
## 1950 1950.000 1950.083 1950.167 1950.250 1950.333 1950.417 1950.500
## 1951 1951.000 1951.083 1951.167 1951.250 1951.333 1951.417 1951.500
## 1952 1952.000 1952.083 1952.167 1952.250 1952.333 1952.417 1952.500
## 1953 1953.000 1953.083 1953.167 1953.250 1953.333 1953.417 1953.500
## 1954 1954.000 1954.083 1954.167 1954.250 1954.333 1954.417 1954.500
## 1955 1955.000 1955.083 1955.167 1955.250 1955.333 1955.417 1955.500
## 1956 1956.000 1956.083 1956.167 1956.250 1956.333 1956.417 1956.500
## 1957 1957.000 1957.083 1957.167 1957.250 1957.333 1957.417 1957.500
## 1958 1958.000 1958.083 1958.167 1958.250 1958.333 1958.417 1958.500
## 1959 1959.000 1959.083 1959.167 1959.250 1959.333 1959.417 1959.500
## 1960 1960.000 1960.083 1960.167 1960.250 1960.333 1960.417 1960.500
##           Aug      Sep      Oct      Nov      Dec
## 1949 1949.583 1949.667 1949.750 1949.833 1949.917
## 1950 1950.583 1950.667 1950.750 1950.833 1950.917
## 1951 1951.583 1951.667 1951.750 1951.833 1951.917
## 1952 1952.583 1952.667 1952.750 1952.833 1952.917
## 1953 1953.583 1953.667 1953.750 1953.833 1953.917
## 1954 1954.583 1954.667 1954.750 1954.833 1954.917
## 1955 1955.583 1955.667 1955.750 1955.833 1955.917
## 1956 1956.583 1956.667 1956.750 1956.833 1956.917
## 1957 1957.583 1957.667 1957.750 1957.833 1957.917
## 1958 1958.583 1958.667 1958.750 1958.833 1958.917
## 1959 1959.583 1959.667 1959.750 1959.833 1959.917
## 1960 1960.583 1960.667 1960.750 1960.833 1960.917
#This will print the cycle across years.
frequency(AirPassengers)
## [1] 12

Step 3:

It is essential to analyze the trends prior to building any kind of time series model. The details we are interested in pertains to any kind of trend, seasonality or random behaviour in the series. what better way to do so than visualize the Time Series.

#This will plot the time series
ts.plot(AirPassengers, xlab="Year", ylab="Number of Passengers", main="Monthly totals of international airline passengers, 1949-1960")

# This will fit in a line
abline(reg=lm(AirPassengers~time(AirPassengers)))

AUTO CORRELATION:

It is a very powerful tool for Time series analysis.Process with auto correlation are more predictable as compared to none.

Total Correlation Chart (also known as Auto – correlation Function / ACF) - auto correlation defined as a function of the time lag. Autocorrelations or lagged correlations are used to assess whether a time series is dependent on its pastFor a time series x of length n we consider the n-1 pairs of observations one time unit apart. The first such pair is (x[2],x[1]), and the next is (x[3],x[2]). Each such pair is of the form (x[t],x[t-1]) where t is the observation index, which we vary from 2 to n in this case. The lag-1 autocorrelation of x can be estimated as the sample correlation of these (x[t], x[t-1]) pairs.

If the above doesnt make sense, luckily, the acf() command provides a shortcut. Applying acf(…, lag.max = 1, plot = FALSE) to a series x automatically calculates the lag-1 autocorrelation.

ACF help us determine what type of series we have, whether it is a White noise, Random walk, Auto regressive or Moving average.

acf(AirPassengers)

Fit the AR model to AirPassengers

For a given time series x we can fit the autoregressive (AR) model using the arima() command and setting order equal to c(1, 0, 0). Note for reference that an AR model is an ARIMA(1, 0, 0) model.

#Fitting the AR Model to the time series
AR <- arima(AirPassengers, order = c(1,0,0))
print(AR)
## 
## Call:
## arima(x = AirPassengers, order = c(1, 0, 0))
## 
## Coefficients:
##          ar1  intercept
##       0.9646   278.4649
## s.e.  0.0214    67.1141
## 
## sigma^2 estimated as 1119:  log likelihood = -711.09,  aic = 1428.18
#plotting the series along with the fitted values
ts.plot(AirPassengers)
AR_fit <- AirPassengers - residuals(AR)
points(AR_fit, type = "l", col = 2, lty = 2)

Forcasting using AR model

The predict() function can be used to make forecasts from an estimated AR model. In the object generated by your predict() command, the $pred value is the forceast, and the $se value is the standard error for the forceast. To make predictions for several periods beyond the last observations, you can use the n.ahead argument in your predict() command. This argument establishes the forecast horizon (h), or the number of periods being forceast. The forecasts are made recursively from 1 to h-steps ahead from the end of the observed time series.

#Using predict() to make a 1-step forecast
predict_AR <- predict(AR)

#Obtaining the 1-step forecast using $pred[1]
predict_AR$pred[1]
## [1] 426.5698
#ALternatively Using predict to make 1-step through 10-step forecasts
predict(AR, n.ahead = 10)
## $pred
##           Jan      Feb      Mar      Apr      May      Jun      Jul
## 1961 426.5698 421.3316 416.2787 411.4045 406.7027 402.1672 397.7921
##           Aug      Sep      Oct
## 1961 393.5717 389.5006 385.5735
## 
## $se
##           Jan      Feb      Mar      Apr      May      Jun      Jul
## 1961 33.44577 46.47055 55.92922 63.47710 69.77093 75.15550 79.84042
##           Aug      Sep      Oct
## 1961 83.96535 87.62943 90.90636
#plotting the AirPassenger series plus the forecast and 95% prediction intervals
ts.plot(AirPassengers, xlim = c(1949, 1961))
AR_forecast <- predict(AR, n.ahead = 10)$pred
AR_forecast_se <- predict(AR, n.ahead = 10)$se
points(AR_forecast, type = "l", col = 2)
points(AR_forecast - 2*AR_forecast_se, type = "l", col = 2, lty = 2)
points(AR_forecast + 2*AR_forecast_se, type = "l", col = 2, lty = 2)

Fit the MA model to AirPassengers

We can fit the simple moving average (MA) model using arima(…, order = c(0, 0, 1)). Note for reference that an MA model is an ARIMA(0, 0, 1) model.

#Fitting the MA model to AirPassengers
MA <- arima(AirPassengers, order = c(0,0,1))
print(MA)
## 
## Call:
## arima(x = AirPassengers, order = c(0, 0, 1))
## 
## Coefficients:
##          ma1  intercept
##       0.9642   280.6464
## s.e.  0.0214    10.5788
## 
## sigma^2 estimated as 4205:  log likelihood = -806.43,  aic = 1618.86
#plotting the series along with the MA fitted values
ts.plot(AirPassengers)
MA_fit <- AirPassengers - resid(MA)
points(MA_fit, type = "l", col = 2, lty = 2)

Forcasting using MA model

#Making a 1-step forecast based on MA
predict_MA <- predict(MA)

#Obtaining the 1-step forecast using $pred[1]
predict_MA$pred[1]
## [1] 425.1049
#Alternately Making a 1-step through 10-step forecast based on MA
predict(MA,n.ahead=10)
## $pred
##           Jan      Feb      Mar      Apr      May      Jun      Jul
## 1961 425.1049 280.6464 280.6464 280.6464 280.6464 280.6464 280.6464
##           Aug      Sep      Oct
## 1961 280.6464 280.6464 280.6464
## 
## $se
##           Jan      Feb      Mar      Apr      May      Jun      Jul
## 1961 64.84895 90.08403 90.08403 90.08403 90.08403 90.08403 90.08403
##           Aug      Sep      Oct
## 1961 90.08403 90.08403 90.08403
#Plotting the AIrPAssenger series plus the forecast and 95% prediction intervals
ts.plot(AirPassengers, xlim = c(1949, 1961))
MA_forecasts <- predict(MA, n.ahead = 10)$pred
MA_forecast_se <- predict(MA, n.ahead = 10)$se
points(MA_forecasts, type = "l", col = 2)
points(MA_forecasts - 2*MA_forecast_se, type = "l", col = 2, lty = 2)
points(MA_forecasts + 2*MA_forecast_se, type = "l", col = 2, lty = 2)

Choosing AR or MA: Exploiting ACF plots

Once we have got the models ready we must answer the important question: Should we choose AR or MA process? Goodness of fit such as an Information criteria is an method to help us make the decision. Specifically Akaike information criterion (AIC) and Bayesian information criterion (BIC) are used for Time series Models. Information Criteria is a more advanced concept but for either measure a lower value indicates a relatively better fitting model.

While the math underlying the AIC and BIC is beyond the scope of this vignettw, for your purposes the main idea is these these indicators penalize models with more estimated parameters, to avoid overfitting, and smaller values are preferred. All factors being equal, a model that produces a lower AIC or BIC than another model is considered a better fit.

# Find correlation between AR_fit and MA_fit
cor(AR_fit, MA_fit)
## [1] 0.954995
# Find AIC of AR
AIC(AR)
## [1] 1428.179
# Find AIC of MA
AIC(MA)
## [1] 1618.863
# Find BIC of AR
BIC(AR)
## [1] 1437.089
# Find BIC of MA
BIC(MA)
## [1] 1627.772

Given the lower value of AIC and BIC in AR model, we should go with that for the time series analysis of AirPassenger data.

Conclusion:

Whilst there are many models applicable to time series data AR and MA are the most common ones. Choosing AR or MA clearly depends upon the type of data we are trying to analyse.