Welcome back to
Practical Time Series Analysis. We've looked in previous
lectures at mixed models, the so-called ARMA models,
Auto-Regressive Moving Average. We've learned to work with
the difference polynomials in order to move between infinite order moving
average, infinite order auto regressive, or order PQ ARMA models. It's time to look at a dataset. When you're done with this lecture,
you should be able to use R to estimate coefficients in an ARMA model, have
a decent understanding of what this is all about, and use a criterion such as
the AKAIKE information criterion or some other related criterion
to select a good model. What is an ARMA model? If we want to understand the state of
the system at time t, we'll build it up as some noise coming in from the environment
plus an auto regressive part. So this would be a weighted average
of past values of the system, plus a moving average part. A weighted average of past influences
on the system through noise. The dataset we would like to use is built
into R, it's called Discoveries, and it looks at scientific inventions. From the middle 1800s to the middle
of 1900s, how many inventions and discoveries did you get in a given year? We'll do the usual sort of plotting. We'll plot the time series sequentially,
and we'll also use stripchart. This is a nice routine to use
when you have discrete data. Rather than a histogram, strip chart will break the data out
by location a little bit better. You can see the time series on the left. And it does bounce around a bit, but
there seems to be some structure here. It doesn't appear to be just noise. If we look aggregated at the sort of
distribution for the number discoveries, you might think about that
as a Poisson distribution. Or if you have a different discreet
distribution, you believe applies and you'd like to fit,
you can chase that down. What we're going to do is try
to build an ARMA model here. In Shumway and Stoffer's book,
Time Series Analysis and Its Applications, there is a table that you'll actually see
in various places on the Internet as well. And it summarizes nicely how we like
to think about our model building. You can obtain the book for free on the Internet as well,
it's a terrific resource. If you have an auto-regressive model,
just pure auto-regressive, then we've seen this before when
we did our theoretical work. The ACF, the auto-correlation function,
will tail off, but the PACF will cut off rather
abruptly after lag p. So that can help us to
determine the order there. Similarly, for a moving average model,
the ACF is going to cut off rather sharply after q, but
the PACF will tail off. Finally, if you have a mix model, both
of these plats will appear to tail off. Now when we do mathematics,
we like have strict inequalities, we like to have results that are clear and
unambiguous. If you have a time series, especially
if it's a little bit on the short side, and especially if it has noise,
these bits of advice, let's say, are harder to apply in
practice than one might think. So, tails off, intuitively we would
believe that tails off means something like exponential decay. Cuts off after lag p, okay, so
that's going to be more abrupt. But if your system is noisy, telling
the difference between something tailing off and cutting off may not
be as easy as you think. So let's see what our plots look like. We'll make the usual calls. The ACF on discoveries and the PACF on discoveries, and there you go. There appear to be two, maybe three
spikes above noise over here on the ACF. That doesn't look like a very
abrupt cutting off to me, seems more like tailing off. But again,
there's some noise in the system. It's a little bit hard to tell. The PACF seems to tail off,
so hard to say here. I am seeing one, two maybe if you're
generous, spikes above noise on the PACF. And it looks like we have
three spikes above noise here. So it's probably a good idea in a case
like this to try several competing models. We'll put in different orders of p and q. We have ways to assess relative quality
of models, the AIC for instance. And we'll do that for a reasonable number of model types or candidate models. So as we explore our ARMA models,
let's let q go from 0 to 3, and p go from 0 to 3 and see what we get. We'll assess quality of
the different models. When I did this, it's a little tedious, the full results are available
in the readings. But you should take five minutes and
do this yourself. You can either write a loop if you're so
inclined, or just make continual calls to ARMA, letting p and q both go from 0 to 3. There just aren't that many of them. When we do that, the AIC,
if you wish to interrogate your model, you can extract the AKAIKE information
criterion with the AIC call. The best one is 440.2 or so, with a 1,0,1 model. We can do very marginally better
with a more complicated model, but these numbers are really very close. Then you have to ask yourself whether
the extra complexity in model is justified by the very
small difference in AIC. When I make a full call to
ARMA with the order 1,0,1, you can see the coefficients
that we obtain and you can see the output that
we get from the ARMA call. Now, there are automatic
routines available. If you have low orders of p and q,
you may just want to do it all yourself. But if you're rushed or lazy, auto.arima
will come to the rescue for you. The only little rub
here is that auto.arima is going to make an approximation. If your time series gets long, then the internal calculations
starts to take a bit of time. So, you can unpack this
little zen co-on down here to tell whether we're going to
use an approximation or, rather, whether auto.arima will use
an approximation on our dataset. Auto.arima lives in the forecast library,
so we'll make a call there. And we didn't see any evidence of trend,
so we're going to let d equal 0. I ran it with approximation=FALSE,
and auto.arima, like say 2,0,0 model,
a second order auto-regressive model. You can see that there are various ways
implied by the output to auto.arima that it will allow the user to make
a choice on various selection criteria. So there's the AIC, the corrected AIC and the Bayesian Information Criterion
are all available to you. You may specify these in the call. And this screen will show
you how it might do that. You can set a flag ic=bic or
aic or the corrected aic. When you send it forward with
the Bayesian information criterion, you return a 1,0,1 model. The AIC returns a 3,0,0, and so on. At this point, we're able to estimate
coefficients in an ARMA model. We can do this with setting calls to
ARMA ourselves with orders p and q. Or if you prefer, you can use you
can use the auto.arima routine. Within auto.arima,
you can specify your selection criterion. So, you get a little flexibility that way. And then you'll use this to try
to build a parsimonious model, to understand your dataset.