ARMA Estimation Example

In this notebook we will fit several \(AR\), \(MA\) and \(ARMA\) models to inflation data. Let us begin by pulling CPI inflation data into R and taking first differences (we will see that the first differences appear to be more stationary than levels).

library(Ecdat)
data(Mishkin)
inflation = Mishkin[,1]
diffInflation = diff(inflation)

To get a sense of the data, let’s plot the time series of inflation and first differences of inflation.

par(mfrow=c(2,1))
plot(inflation, ylab="Inflation rate", main="(a)")
abline(h=0)

plot(diffInflation, ylab="Change in inflation rate", main="(b)")
abline(h=0)

Inflation levels do not appear to have a stationary mean, however the first differences do. The autocorrelations of first differences are depicted below.

acf(as.numeric(diffInflation))

The ACF plot suggests that an \(MA(1)\) model might be appropriate for this data. We will begin, however, by fitting an \(AR(1)\) and examining the residuals.

ar1Fit = arima(diffInflation, order=c(1,0,0));
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
plot(ar1Fit$resid)
abline(h=0)

acf(as.numeric(ar1Fit$resid))
qqnorm(ar1Fit$resid, datax=TRUE)

qqline(ar1Fit$resid, datax=TRUE)

The time series and ACF plots of the residuals look decent, although the negative autocorrelations at lags 2 and 3 might suggest that we still haven’t appropriately modeled the time dynamics of the data. The next code block fits an \(AR(p)\) model for \(p = 1,2,\ldots,20\) and plots both the AIC and BIC values for each fit.

#install.packages("forecast"); # run this if 'forecast' not installed
library(forecast)
pMax = 20
aicAR = bicAR = rep(0, length=pMax)
for(p in 1:pMax){
    out = arima(diffInflation, order=c(p,0,0))
    aicAR[p] = out$aic
    bicAR[p] = out$aic - 2*(2+p) + log(length(diffInflation))*(2+p)
}
yMin = min(min(aicAR), min(bicAR))
yMax = max(max(aicAR), max(bicAR))
plot(aicAR, xlab="p", ylab="Criterion", pch=8, ylim=c(yMin,yMax),
    main="AR Fits to Changes in Inflation Rate")
points(bicAR)

legend("topright", legend=c("AIC", "BIC"), pch=c(8,1))

The plot shows that BIC is minimized at \(p = 6\) while AIC is minimized at \(p = 16\). We can also use the auto.arima function to automatically fit and select the best \(AR\) model using AIC or BIC. The following code does this and then examines the residuals for the best fit model.

out = auto.arima(diffInflation, max.p=10, max.q=0, ic="bic")
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
plot(out$resid)
abline(h=0)

acf(out$resid)
qqnorm(out$resid, datax=TRUE)

qqline(out$resid, datax=TRUE)

The ACF of the residuals looks good. Let’s now repeat the analysis above, fitting an \(MA(1)\) and then plotting the AIC and BIC for a sequence of \(MA(q)\) models, for \(q = 1,2,\ldots,10\).

ma1Fit = arima(diffInflation, order=c(0,0,1))
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
plot(ma1Fit$resid)
abline(h=0)

acf(ma1Fit$resid)
qqnorm(ma1Fit$resid, datax=TRUE)

qqline(ma1Fit$resid, datax=TRUE)

qMax = 10
aicMA = bicMA = rep(0, length=qMax)
for(q in 1:qMax){
    out = arima(diffInflation, order=c(0,0,q))
    aicMA[q] = out$aic
    bicMA[q] = out$aic - 2*(2+q) + log(length(diffInflation))*(2+q)
}
yMin = min(min(aicMA), min(bicMA))
yMax = max(max(aicMA), max(bicMA))
plot(aicMA, xlab="q", ylab="Criterion", pch=8, ylim=c(yMin,yMax),
    main="MA Fits to Changes in Inflation Rate")
points(bicMA)

legend("topright", legend=c("AIC", "BIC"), pch=c(8,1))

The BIC is minimized at \(q=3\) and the AIC is minimized at \(q=2\). We can now use the auto.arima function to automatically select the best \(MA\) model and to examine the corresponding residuals.

out = auto.arima(diffInflation, max.p=0, max.q=10, ic="bic")
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
plot(out$resid)
abline(h=0)

acf(out$resid)
qqnorm(out$resid, datax=TRUE)

qqline(out$resid, datax=TRUE)

Finally, we use the auto.arima function to fit a sequence of \(ARMA(p,q)\) models for \(p=\{1,2,3,4,5\}\) and \(q=\{1,2,3,4,5\}\) and examine the residuals of the best fit model.

out = auto.arima(diffInflation, max.p=5, max.q=5, ic="bic")
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
plot(out$resid)
abline(h=0)

acf(out$resid)
qqnorm(out$resid, datax=TRUE)

qqline(out$resid, datax=TRUE)