Transforming imaginary and complex numbers

Difference equations are the topic of interest in Chapter 1 of Hamilton's Time Series Analysis. The chapter includes a section that outlines the general solution of pth-order difference equations with distinct eigenvalues. On page 14, Hamilton says that:

interesting possibilities arise if some of the eigenvalues are complex

and then a series of equations follow this statement. To guide the reader through the equations, footnotes are provided at the bottom of the page that refer the reader to the mathematical appendix at the back of the book.

It can be a little bit difficult to figure out what's going on in the equations on pages 14 to 16 of the book, so I thought it would be a good idea to write up some notes in order to make the material more understandable.

The equations involve transforming imaginary and complex numbers to trigonometric functions. It's straightforward stuff really, but since no reference is made to the names of the formulas being used, it may leave you scratching your head as to what's going on.

So, our aim here is to understand how to transform imaginary and complex numbers to trigonometric functions. To do this, we're going to need to familiarize ourselves with the following:

1. Expressing Cartesian coordinates in polar coordinate form; 
2. Euler's formula; 
3. De Moivre's Theorem.

Let's begin with expressing Cartesian coordinates in polar coordinate form.

$a$ and $b$ in the diagram below are called Cartesian coordinates. These can be expressed in terms of $\theta$ and $R = \sqrt{a^{2}+b^{2}}$ which are called polar coordinates. To do this, we use the following formulas: \[a = R~cos~\theta,~ ~~~~~ b = R~sin~\theta,~ ~~~~~ R>0.\] Applying these formulas to the conjugate complex number $(a \pm b i)$, we have: \[a \pm b i = R~cos~\theta \pm i~R~sin~\theta = R(cos~\theta \pm i~sin~\theta) \] and this explains Equation [1.2.36] in Hamilton.

Second on the list is Euler's formula. This tells us that: \[ e^{\pm i \theta} = cos~\theta \pm i~sin~\theta.\] Thus we can express $(a \pm bi) = R(cos~\theta \pm i~sin~\theta)$ as \[a \pm bi = R(cos~\theta \pm i~sin~\theta) = R e^{\pm i \theta}.\] The last item on the list is De Moivre's theorem, but first let's see how a conjugate complex number is raised to the jth power. It is done as follows: \[(a \pm bi)^{j} = (Re^{\pm i \theta})^{j} = R^{j} e^{\pm i j \theta}.\] Or, using De Moivre's theorem we have: \[(a \pm bi)^{j} = R^{j}(cos~j~\theta \pm i~sin~j~\theta).\] and now we've arrived at Equation [1.2.37] in Hamilton.

Although not explicitly mentioned in the textbook, it's Euler's formula and De Moivre's theorem that lead to both Equation [1.2.37] and Equation [1.2.38] in Chapter 1 of Hamilton.

The rest of the section on the general solution of pth-order difference equations with distinct eigenvalues involves no more than basic substitution and a little bit of algebra. Despite the fact that the notation looks a bit "heavy", it's easy to follow once you know to use Euler's formula and De Moivre's theorem.

Impulse Response Functions

I bought a copy of Ed Leamer's book Macroeconomic Patterns and Stories after hearing him speaking about it on EconTalk. Actually, Professor Leamer is probably my favourite guest on EconTalk, so purchasing his book was a no-brainer for me. If you've heard of him before it's probably due to his famous paper Let's Take the Con Out of Econometrics. He is also well known for advocating "extreme bounds analysis", which, in basic terms, involves finding a possible range of values for an estimate and investigating the fragility of the model by seeing if the value of the estimate is sensitive to small changes in model specification.

The example of EBA given in Let's Take the Con Out of Econometrics is to do with the effect of capital punishment on the murder rate. Using a number of regressions, with different, but somewhat similar specifications, Leamer shows that at one end of the spectrum, the estimate is negative, but at the other end, the estimate is positive!

In his book, Macroeconomic Patterns and Stories, Professor Leamer talks about four key macro variables in detail; namely, Growth, Unemployment, Inflation, and Interest Rates. In this post, I want to focus on these four variables and in particular to the topic of Chapter 6 of the book which is Extrapolative Forecasting.

In Chapter 6, a table is shown and it contains regression results. There are four separate regressions and each of them are autoregressive. Growth is modeled as an AR(2) process, so too is the unemployment rate and 3-month treasuries, while inflation is modeled as an AR(4) process. Admittedly, Professor Leamer says that a more careful data analysis could be done on the data, but the results from these regressions still permit us to get an insight into the dynamic nature of the four variables.

In R, I've estimated how each of the four variables respond to a one unit shock. I've plotted the results in the Figure 1 below. A similar graph appears on page x in Leamer's book and this gives me confidence that the code I wrote in R is correct.

Figure 1

To explain the dynamic behaviour of each of the four variables in response to a shock, Leamer introduces two concepts to the reader: persistence and momentum.

R Time Series Issue: Residuals (Revisited)

In the previous post, I investigated the arima() function and its method of calculating residuals. The issue that I raised was that the first residual associated with arima(x, c(1,0,0)) happens to be non-zero when it should probably either be treated as zero or classified as n.a. (I prefer the latter) when doing time series analysis.

Last time, I specified only two arguments in the arima() function; namely, x and order. For more information type ?arima in R.

In this post, I'd like to point out some more details about calculating residuals in R via arima() and to do so, I'll talk about a third argument of arima() - the method argument - which is used to specify which estimation method is to be used to estimate the parameters of the model. I briefly mentioned last time that maximum-likelihood and non-linear least squares were two possible estimation methods for ARIMA models.

The method argument has three options: "CSS-ML", "ML", and "CSS". The default is set to "CSS-ML", so that was the method used in the previous post. The help page for the arima() function states the following about the method argument:

Fitting method: maximum likelihood or minimize conditional sum-of-squares. The default (unless there are missing values) is to use conditional-sum-of-squares to find starting values, then maximum likelihood. 

I cannot confirm it yet, but I'm almost certain that "CSS-ML" (the default) refers to Marquardt's compromise, "ML" to maximum likelihood, and "CSS" to what I (and the documentation that I've read) call non-linear least squares.

To begin my investigation, I estimated an AR(1) model for the same data in the previous post using, in turn, each of the three estimation methods. That is, I estimated these three models in R:

1. arima(x, c(1,0,0), method = c("CSS-ML")),
2. arima(x, c(1,0,0), method = c("ML")),
3. arima(x, c(1,0,0), method = c("CSS")).

Now, here is the problem...

Upon inspection of the residuals associated with each model, I found that, in two of the three cases ("CSS-ML", "ML") the first residual is non-zero, while in the third case ("CSS"), the first residual is precisely equal to zero.

I've tabulated the residuals for each of the three models in Table 1 below...

If the analyst accepts these first residuals as appropriate then they can fire away and use x$\$$residuals as they please. On the other hand, if the analyst is of the opinion that the first residual should be omitted from their diagnostic checking procedures, they can do this in R by specifying x$residuals[-1]. This command will return all but the first residual as a vector object.

R Time Series Issue: Residuals

Professor David S. Stoffer at the Department of Statistics at the University of Pittsburgh has documented some issues with using R for time-series analysis. In this post, I would like to make an addition to the issues pointed out by Professor Stoffer.

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.

Table 1 \begin{array}{|c|c|} \hline t & z_{t} \\ \hline 0 & n.a. \\ 1 & 4.4 \\ 2 & 5.8 \\ 3 & 6.7 \\ 4 & 7.1 \\ 5 & 5.7 \\ \vdots & \vdots \\ 59 & 13.3 \\ 60 & 6.2 \\ \hline \end{array} Now suppose that given this data, we use the Box-Jenkins methodology and happen to identify an AR(1) process as a tentative model. That is, we investigated the acutocorrelation function and partial autocorrelation for the series and found that the autocorrelation decayed rapidly and there was a single significant spike in the partial autocorrelation function.

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:

Table 2 \begin{array}{|c|c|} \hline t & \hat{a}_{t} \\ \hline 1 & -1.2037423 \\ 2 & 0.8749301 \\ 3 & 0.8225200 \\ 4 & 0.6102564 \\ \vdots & \vdots \\ 60 & -4.7796768 \\ \hline \end{array} Where is this extra residual that R calculates, namely $\hat{a}_{1} = -1.2037423$, coming from? Recall that we do not calculate $\hat{z}_{1}$ and subsequently $\hat{a}_{1}$ because these calculations involve $z_{0}$ which we do not know. See Table 1, it's n.a.!

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.

How to install the tseries package

I recently installed the tseries package because I wanted to make use of its bds.test() function. This function allows one to carry out the BDS test. It's one of the most popular tests for testing the iid hypothesis and it can also be used to test for non-linearity.

When I tried to install the tseries package in R using install.packages(tseries), I received an error message saying that the package could not be found or something like that there was no such package. Installing packages in R is usually straight forward enough, so I was slightly puzzled on this occasion.

I figured my best bet would be to go directly to the CRAN website and download the tseries source package tarball file then install it using the Terminal in Ubuntu. After trying to install the package from the terminal, I received a second error message. This time I was informed that I'd need to first download and install another package called quadprog.

So, I downloaded and installed the quadprog package and then finally this enabled me to successfully install the tseries package. Some others have encountered the same problem when trying to install the tseries package, so in case you do to,  you may find it useful to know the commands that I used to get it all sorted out.

If you're not used to using the terminal, you may like to know that the wget command is a utility for non-interactive download of files from the Web - you simply type wget followed by the url of the file to be downloaded. The sudo command basically allows you to execute commands as a superuser or administrator on your computer. If you need to, you can read about R CMD INSTALL and how to install R packages here.