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.