The issue I'd like to bring to the fore is on the subject of residuals and how they are calculated after using the arima() function. Let's go through a numerical example - I suggest grabbing a pencil and paper - with an application in R:
Suppose we have a realization $ z_{t} = (z_{1}, z_{2}, \ldots, z_{t-1}, z_{t})$ where $t = 60$. The data is shown in Table 1 below, although I have omitted a number of data points because we won't be needing them here (don't worry, you can download all of the data, and an R script, at the bottom of the post in order to replicate all of the results presented here). Notice that I have included a row in the table showing that we do not have an observation for $z_{t}$ where $t = 0$. This will be important later on and will, in fact, show us the crux of the problem that I'm raising.
In backshift notation an AR(1) can be written as \[(1 - \phi_{1}B) \tilde{z}_{t} = a_{t} \] where $\tilde{z}_{t} = z_{t} - \mu$ is just the original time-series expressed in deviations from the mean, $a_{t}$ are random shocks, B is the backshift operator, and $(1 - \phi_{1}B)$ is an AR operator which can be written in a generalized AR(p) form as: $(1 - \phi_{1}B - \phi_{2}B^{2} - \cdots - \phi_{p}B^{p})$.
Rewriting the model in common algebraic form, because this will help with explanations later on, results in the following series of equations: \begin{aligned} (1 - \phi_{1}B) \tilde{z}_{t} &= a_{t}\\ (1 - \phi_{1}B) (z_{t} - \mu) &= a_{t}\\ z_{t} - \mu - \phi_{1}B z_{t} + \phi_{1}B \mu &= a_{t}\\ z_{t} - \mu - \phi_{1} z_{t-1} + \phi_{1} \mu &= a_{t}\\ z_{t} &= \mu + \phi_{1} z_{t-1} - \phi_{1} \mu + a_{t}\\ z_{t} &= \mu(1 - \phi_{1}) + \phi_{1} z_{t-1} + a_{t}\\ z_{t} &= C + \phi_{1} z_{t-1} + a_{t} \end{aligned} where $\mu(1 - \phi_{1}) = C$ is the constant term. (Tip: Chapter 2 of Hamilton's Time Series Analysis gives an excellent round-up of the algebra of the backshift (or lag) operator. Note, however, that I've used the letter B where Hamilton uses the letter L. Additionally, a concise list of properties of the lag operator can be found in Chapter 8 (p.168-9) of Harrison & Waldron's Mathematics for Economics and Finance.)
At the second stage of the Box-Jenkins methodology - the estimation stage - we fit the model and find estimates for the parameters in our model. In this case, we need to find estimates for $\phi_{1}$ and the constant, $C$. These estimates of the true underlying values are denoted $\hat{\phi}_{1}$ and $\hat{C}$, respectively. ARIMA estimation is usually carried out using nonlinear least squares, however, there is a trade off between using least squares and maximum likelihood which will not discussed further since it is outside the scope of this post.
In R, we can fit an AR(1) model using arima(x, order=c(1,0,0)) and in doing so using the data from Table 1, we get the following results: \[ \hat{\phi}_{1} = 0.6803 \text{ and } \hat{C} = 1.931723.\] Actually, R reports $\hat{\mu}$ as the intercept, so it was necessary to calculate $\hat{C}$ "by hand". That is, the arima() output tells us that the intercept is $6.0423$, but when it refers to the intercept it means $\hat{\mu}$ rather than $\hat{C}$. This issue has already been documented by Professor Stoffer, but to make it clear, R reports to us the following information: \begin{aligned} z_{t} &= \mu(1 - \phi_{1}) + \phi_{1} z_{t-1} + a_{t}\\ z_{t} &= 6.0423(1 - 0.6803) + 0.6803 z_{t-1} + a_{t}\\ z_{t} &= 1.931723 + 0.6803 z_{t-1} + a_{t} \end{aligned} We are NOT told, as one may suspect, that: \begin{aligned} z_{t} &= C + \phi_{1} z_{t-1} + a_{t}\\ z_{t} &= 6.0423 + 0.6803 z_{t-1} + a_{t} \end{aligned}
To calculate predicated values, denoted $\hat{z}_{t}$ (that's a "z hat" and not a "z tilde" which I already said represents deviations from the mean), we simply plug values into the following equation. Fortunately, this is done automatically in R, but here are examples using our estimates and the data shown in Table 1 above: \begin{aligned} \hat{z}_{t} &= \hat{C} + \hat{\phi}_{1} z_{t-1}\\ \hat{z}_{2} &= 1.931723 + 0.6803 (4.4) = 4.925043 \\ \hat{z}_{3} &= 1.931723 + 0.6803 (5.8) = 5.877463 \\ \hat{z}_{4} &= 1.931723 + 0.6803 (6.7) = 6.484733 \\ \hat{z}_{60} &= 1.931723 + 0.6803 (13.3) = 10.979713 \end{aligned} Notice that we do not calculate $\hat{z}_{1}$ because in that case we have $\hat{z}_{1} = \hat{C} + \hat{\phi}_{1} z_{0}$ and data is not available for $t=0$.
To calculate the residuals, denoted $\hat{a}_{t}$, for our AR(1) model we subtract the predicted values from the observed values. Again, this is done automatically in R, but here again are some examples: \begin{aligned} \hat{a}_{t} &= z_{t} - \hat{z}_{t}\\ \hat{a}_{2} &= 5.8 - 4.925043 = 0.874957\\ \hat{a}_{3} &= 6.7 - 5.877463 = 0.827537\\ \hat{a}_{4} &= 7.1 - 6.484733 = 0.615267\\ \hat{a}_{60} &= 6.2 - 10.979713 = -4.779713 \end{aligned} It has been shown here that estimating an AR(1) model using a realization with $t=60$ gives rise to 59 residuals. Now here is the problem...
Let's take a look at the residuals that are calculated in R using arima(). To do this, we can assign the model arima(x, c(1,0,0)) using our data to an object called x and then check out the attributes of object x; type in R attributes(x). We find that there is an attribute name called "residuals". These are the residuals of the estimated model. Viewing the residuals using x$residuals we find that 60 residuals have been calculated, not 59 as we should expect given the calculations performed above! R reports to us the residuals shown in Table 2:
To decipher what R appears be doing, consider the following: \begin{aligned} \hat{z}_{t} &= \hat{C} + \hat{\phi}_{1} z_{t-1}\\ \hat{z}_{1} &= \hat{C} + \hat{\phi}_{1} z_{0} \end{aligned} If we treat $z_{0}$ as zero we end up with: \[\hat{z}_{1} = \hat{C} = 1.931723\] Now we can calculate a residual for $t=1$ in this way: \begin{aligned} \hat{a}_{t} &= z_{t} - \hat{z}_{t}\\ \hat{a}_{1} &= 4.4 - 1.931723 = 2.468277 \end{aligned} But this value certainly does not correspond to the value of $\hat{a}_{1}$ that R calculates; namely $\hat{a}_{1} = -1.2037423$.
OK, let's go back and again consider the following: \[\hat{z}_{1} = \hat{C} + \hat{\phi}_{1} z_{0}\] but let's rewrite it and work with the equation in this form: \[\hat{z}_{1} = \mu(1 - \phi_{1}) + \phi_{1} z_{0}.\] If we again treat $z_{0}$ as zero we end up with: \begin{aligned} \hat{z}_{1} &= \mu(1 - \phi_{1}) \\ \hat{z}_{1} &= 6.0423(1 - 1.931723) = -5.629749 \end{aligned} Calculating the residual using this predicted value, we get \[\hat{a}_{1} = 4.4 - 5.629749 = -1.229749\] which is close to R's calculation of $\hat{a}_{1}$, namely $\hat{a}_{1} = −1.2037423$.
I'm thus led to believe that this is what's going on when R reports this extra residual for $t=1$. I have tried a few other procedures, but they point in the wrong direction and end up way off the mark.
The question now to be raised is that when using R for time-series analysis, at the diagnostic stage of the Box-Jenkins methodology, should the analyst exclude this extra (non-zero) residual or keep it when performing diagnostic checks? Furthermore, ask yourself what happens to the residuals in R when an AR(2), say, with 60 observations is estimated.
Note: expect some small discrepancies between my calculated figures and R's calculated figures, this is down to rounding and shouldn't be much of an issue here.
No comments:
Post a Comment