Saturday, 25 May 2013
CRAN Task Views
In the previous post, I provided a list of R packages that I find useful for time-series analysis, economic and financial research. There are actually expert summaries of R packages by subject on CRAN. This is a really useful resource and I'd highly recommend checking it out.
Thursday, 23 May 2013
R Packages for Econometrics & Time-Series Analysis
During my masters I learned how to use STATA. The STATA documentation is very good, particularly the [TS] manual. I like using STATA, but over time, I came across some tasks that I could not readily perform in STATA. Often, I wanted to make a unusual calculation or estimate a particular model, but I found that these tasks would be difficult or simply impossible to perform in STATA.
Essentially, STATA is more or less a "canned packaged" statistical tool and it has its limitations. There is, of course, the programming language version of STATA called MATA, but at the time I was fairly uninformed about its capabilities. Having since programmed using MATA, I'd still say that the data analysis and modelling tasks that interest me would be difficult to perform in MATA - even for an advanced user.
For the reasons above (and more), I started using the R programming language. R is not a "canned packaged" statistical software; it is a programming language. It is used for statistical analysis, but it has many other uses and different people use it for different reasons. Some people use it for econometrics, others for biometrics, mathematical modelling, statistics, and so on. On the whole, the limitations of R basically come down to the capabilities of the user; if you can program then the sky is the limit.
My general approach to using R and modelling data is that I like to write the code myself. That is, I try not to rely on add-in's or "packages" that have been written by other users of the language. In this way, I know exactly all of the formulas being used in my analysis since I program them myself. This gives me great confidence to be able to stand over the figures and results of my analysis. It would be difficult to ask me a question about where a particular figure is coming from or how it has been made that I could not answer with confidence.
Having said that, there are some really smart people out there writing R packages; I'm talking about leaders of the field in time-series analysis, statistics, etc. So, given that such people's work ought to be deemed reliable, an R user should probably use their packages if they are well documented, tested and trusted.
In other words, writing your own code some of the time is great, but writing it ALL of the time is probably an inefficient way to allocate your time. For instance, if an algorithm has been translated into R code, say, for estimating ARIMA or VAR models then you may as well take advantage of it - not without having read the source documentation, though, of course! Using a black-box for data analysis or forecasting could spell trouble for numerous reasons!
The problem is then, since there are bucket-loads of R packages which ones should you use? Well, here I'd like to suggest a few R packages that I have come across and found useful for time-series analysis and econometrics. In the next few posts I will provide some examples of what they can do.
In no particular order, they are as follows:
Essentially, STATA is more or less a "canned packaged" statistical tool and it has its limitations. There is, of course, the programming language version of STATA called MATA, but at the time I was fairly uninformed about its capabilities. Having since programmed using MATA, I'd still say that the data analysis and modelling tasks that interest me would be difficult to perform in MATA - even for an advanced user.
For the reasons above (and more), I started using the R programming language. R is not a "canned packaged" statistical software; it is a programming language. It is used for statistical analysis, but it has many other uses and different people use it for different reasons. Some people use it for econometrics, others for biometrics, mathematical modelling, statistics, and so on. On the whole, the limitations of R basically come down to the capabilities of the user; if you can program then the sky is the limit.
My general approach to using R and modelling data is that I like to write the code myself. That is, I try not to rely on add-in's or "packages" that have been written by other users of the language. In this way, I know exactly all of the formulas being used in my analysis since I program them myself. This gives me great confidence to be able to stand over the figures and results of my analysis. It would be difficult to ask me a question about where a particular figure is coming from or how it has been made that I could not answer with confidence.
Having said that, there are some really smart people out there writing R packages; I'm talking about leaders of the field in time-series analysis, statistics, etc. So, given that such people's work ought to be deemed reliable, an R user should probably use their packages if they are well documented, tested and trusted.
In other words, writing your own code some of the time is great, but writing it ALL of the time is probably an inefficient way to allocate your time. For instance, if an algorithm has been translated into R code, say, for estimating ARIMA or VAR models then you may as well take advantage of it - not without having read the source documentation, though, of course! Using a black-box for data analysis or forecasting could spell trouble for numerous reasons!
The problem is then, since there are bucket-loads of R packages which ones should you use? Well, here I'd like to suggest a few R packages that I have come across and found useful for time-series analysis and econometrics. In the next few posts I will provide some examples of what they can do.
In no particular order, they are as follows:
- Ecdat: Data sets for econometrics
- pear: Package for Periodic Autoregression Analysis
- vars: VAR Modelling
- strucchange: Testing, Monitoring, and Dating Structural Changes
- dse: Dynamic Systems Estimation (time series package)
- KFAS: Kalman Filter and Smoother for Exponential Family State Space Models
- dlm: Bayesian and Likelihood Analysis of Dynamic Linear Models
- tsDyn: Nonlinear time series models with regime switching
- ccgarch: Conditional Correlation GARCH models
- rugarch: Univariate GARCH models
- dlnm: Distributed Lag Non-linear Models
- urca: Unit root and cointegration tests for time series data
- apt: Asymmetric Price Transmission (apt)
- quantmod: Quantitative Financial Modelling Framework
- FinancialInstrument: Financial Instrument Model Infrastructure for R
- TTR: Technical Trading Rules
- RandomFields: Simulation and Analysis of Random Fields
Saturday, 4 May 2013
Non-linearity in the Irish Stock Market
In this post, I'm going to go through an application of the BDS test using the R programming language. To be specific, I will test for non-linearity in the Irish Stock Market.
The application will involve daily data for the ISEQ index. The data spans from 1999-10-04 to 2013-05-03. There are a total of 3444 observations. Recall from my previous post that for the BDS test to have a reliable performance, the sample size should be greater than 500 observations. The data was obtained from Yahoo! Finance and you can download a .txt file of the data here.
The method that I am going to use involves taking differences of the logs of the data and then applying the BDS test to this transformed series. Note, however, that there is another method of using the BDS test to test for non-linearity which is to fit an ARIMA model (or say a GARCH model) and then to apply the BDS test to the residuals of the model. For further examples of the method that I am going to use in this post, see Chu, P.K.K. (2001) or these course notes by Eric Zivot from the University of Washington.
Let's begin by plotting the data to see what we're working with. The figure below shows a plot of the raw data.
The data is shown in logged form in the next figure. Generally speaking, taking logs has the effect of stabilizing the variance.
Now let's test the data expressed in differences of the logs for non-linearity. We can do this in R by using the bds.test() function from the tseries package. The hypothesis of the test is as follows:
These results suggest that the time-series data for the ISEQ index is non-linearly dependent. This evidence could be used to justify exploring the possibility of building a non-linear model for the series. Keep in mind that we'd normally like to build simpler linear models and not necessarily complicate matters by building non-linear models, but based on the evidence presented here, perhaps a non-linear approach to modelling is appropriate for this particular series.
Here is my R script for this analysis. You can use it to replicate what I've done in this post. I've also included in the script the code used for creating the LaTeX tables.
The application will involve daily data for the ISEQ index. The data spans from 1999-10-04 to 2013-05-03. There are a total of 3444 observations. Recall from my previous post that for the BDS test to have a reliable performance, the sample size should be greater than 500 observations. The data was obtained from Yahoo! Finance and you can download a .txt file of the data here.
The method that I am going to use involves taking differences of the logs of the data and then applying the BDS test to this transformed series. Note, however, that there is another method of using the BDS test to test for non-linearity which is to fit an ARIMA model (or say a GARCH model) and then to apply the BDS test to the residuals of the model. For further examples of the method that I am going to use in this post, see Chu, P.K.K. (2001) or these course notes by Eric Zivot from the University of Washington.
Let's begin by plotting the data to see what we're working with. The figure below shows a plot of the raw data.
The data is shown in logged form in the next figure. Generally speaking, taking logs has the effect of stabilizing the variance.
And here is the data in differences of the logs (which can be interpreted as daily percent change):
- H$_{0}$ : the series is i.i.d.
- H$_{1}$ : the series is not i.i.d. and this implies that non-linearity is present.
BDS Test Statistics
\begin{array}{rrrrr}
\hline
m & \epsilon = 0.5 & \epsilon = 1.0 & \epsilon = 1.5 & \epsilon = 2.0 \\
\hline
2 & 13.81 & 15.47 & 15.59 & 14.15 \\
3 & 19.99 & 21.41 & 21.49 & 20.08 \\
4 & 25.16 & 25.43 & 24.62 & 22.88 \\
5 & 30.00 & 28.94 & 27.12 & 25.07 \\
\hline
\end{array}
p-values
\begin{array}{rrrrr}
\hline
m & \epsilon = 0.5 & \epsilon = 1.0 & \epsilon = 1.5 & \epsilon = 2.0 \\
\hline
2 & 0.00 & 0.00 & 0.00 & 0.00 \\
3 & 0.00 & 0.00 & 0.00 & 0.00 \\
4 & 0.00 & 0.00 & 0.00 & 0.00 \\
5 & 0.00 & 0.00 & 0.00 & 0.00 \\
\hline
\end{array}
The tables above show that all of the test statistics are greater than 1.96, so we should reject the null hypothesis of i.i.d. The p-values provide further evidence for rejecting the null.
These results suggest that the time-series data for the ISEQ index is non-linearly dependent. This evidence could be used to justify exploring the possibility of building a non-linear model for the series. Keep in mind that we'd normally like to build simpler linear models and not necessarily complicate matters by building non-linear models, but based on the evidence presented here, perhaps a non-linear approach to modelling is appropriate for this particular series.
Here is my R script for this analysis. You can use it to replicate what I've done in this post. I've also included in the script the code used for creating the LaTeX tables.
Friday, 3 May 2013
BDS Test for Non-Linearity
A widely used test for non-linearity is the BDS test due to Brock et al. (1987), but later published as Brock et al. (1996). Strictly speaking, the BDS test is a test of independence (a test of the i.i.d hypothesis), but it can be used to test for non-linearity and also as a test of model misspecification.
In this post I'd like to give a brief overview of the BDS test and at the end I'll provide an application of carrying out the BDS test in R. To introduce the BDS test, let me begin by saying that the method of the BDS test has its roots in the concept of spatial correlation from chaos theory.
Let's assume that $\{y_{t}\}$ is a strictly stationary process. Now consider the following vector called the m-history of the process: \[{\bf y}_{t}^{m} = (y_{t}, y_{t-1}, \ldots, y_{t+m-1})\] The construction of the m-history means that the univariate series $y_{t}$ has been embedded into a higher dimensional space. The purpose of this operation is to be able to examine the spatial correlation of the time-series. Notice that the vector, ${\bf y}_{t}^{m}$, consists of overlapping series. This becomes clear if we write out the following equations: \begin{aligned} y_{1}^{m} &= (y_{1}, y_{2}, \ldots, y_{m})\\ y_{2}^{m} &= (y_{2}, y_{3}, \ldots, y_{m+1})\\ & ~\vdots\\ y_{T-m+1}^{m} &= (y_{T-m}, y_{T-m+1}, \ldots, y_{T}) \end{aligned} To demonstrate the concept of embedding (which is due to Dutch mathematician Floris Takens), here is an example from Cromwell et al. (1994). Suppose we have a time-series, $y_{t}$, where $t = 1, 2, \ldots, 5$ and $m = 3$. Embedding this series into higher dimensions results in the following sequence of overlapping series: \begin{aligned} y_{1}^{3} &= (y_{1}, y_{2}, y_{3})\\ y_{2}^{3} &= (y_{2}, y_{3}, y_{4})\\ y_{3}^{3} &= (y_{3}, y_{4}, y_{5}) \end{aligned} In this example, the embedding operation creates three new three dimensional vectors: $y_{1}^{3}, y_{2}^{3}, y_{3}^{3}$. A requirement is that the vectors have the same length. Once the series has been embedded, we need an instrument to allows us to measure (a metric).
The correlation integral, introduced by Grassberger and Procaccia (1983), is a measure of spatial correlation. It measures the distance between the points in the m-dimensional space constructed by the embedding procedure; in the example above, it measures the distance between points in three dimensions. The correlation integral can also be interpreted as a measure of the frequency with which temporal patterns are repeated in the data.
The correlation integral at embedding dimension m is defined as: \[ C_{m, T} (\epsilon) = T_{m} (T_{m} - 1)^{-1} \sum_{t \neq s} I[y_{t}^{m}, y_{s}^{m}, \epsilon ] \] where $T_{m} = T - m + 1$, $t$ and $s$ both range from $1$ to $T - m + 1$ in the summation, m is the embedding dimension, and $\epsilon$ is the threshold value. $I(\cdot)$ is called the indicator function and is defined as: \begin{aligned} I[y_{t}^{m}, y_{s}^{m}, \epsilon ] &= 1 ~~~~~\text{if}~ || y_{t}^{m} - y_{s}^{m} || \leq \epsilon \\ &= 0 ~~~~~\text{otherwise} \end{aligned} The term, $|| y_{t}^{m} - y_{s}^{m} ||$, is a metric for the distance between two pairs. An example using the aforementioned hypothetical data set would be that the distance between $y_{1}^{3} = [y_{1}, y_{2}, y_{3}]$ and $y_{3}^{3} = [y_{3}, y_{4}, y_{5}]$ is given by: \[y_{1}^{3} - y_{3}^{3} = [ y_{1} - y_{3}, y_{2} - y_{4}, y_{3} - y_{5} ] \] The next step of the test involves choosing the maximum element from this resulting vector and then comparing it to the threshold value, denoted, $\epsilon$. If the maximum value is greater than $\epsilon$ then the pair is counted. Otherwise, the pair is not counted. The counting continues as the values for $t$ and $s$ are varied. Then, the number of counted pairs is divided by the total number of possible pairs. This fraction of pairs that are within a distance of $\epsilon$ of each other is the computed value for the correlation integral. The extreme cases are that all pairs are chosen, $C_{m, T} (\epsilon) = 1$, and that no pairs are chosen, $C_{m, T} (\epsilon) = 0$.
The BDS statistic is: \[ BDS = \sqrt{T} \frac{C_{m,T}(\epsilon) - C_{1,T}^{m}(\epsilon)}{\sigma_{m}(\epsilon)}\] where the variance, $\sigma_{m}^{2}(\epsilon)$, has the following complicated definition: \[ \sigma_{m}^{2}(\epsilon) = 4\{K(\epsilon)^{m} + 2 [\sum_{i=1}^{m-1} K(\epsilon)^{m-i} C(\epsilon)^{2i} + (m-1)^{2} C(\epsilon)^{2m} - m^{2} K (\epsilon) C(\epsilon)^{2m-2}]\} \] The indicator functions $C(\epsilon) = E\{I[y_{i}, y_{j}, \epsilon ]\}$ and $K(\epsilon) = E\{I[y_{i}, y_{j}, \epsilon ] I[y_{j}, y_{k}, \epsilon ]\}$ are defined in a similar way. For a complete definition, see Brock et al. (1996).
The null hypothesis of the BDS test is $H_{0}: y_{t}$ is i.i.d. The test is a two tailed test, so we should reject the null hypothesis of i.i.d at the 5% significance level when $|BDS| > 1.96$.
Notice that the statistic is a function of two parameters; the embedding dimension, m, and the threshold value, $\epsilon$. Values for these two parameters must be chosen, so this raises the question of how to choose values for m and $\epsilon$.
Various rules of thumb have been used to choose a value or values for $\epsilon$ and a short list of references can be found be referring to the concluding section of Brock (1995). The default values using the bds.test() function in R are $\epsilon$ = (0.5, 1.0, 1.5, 2.0) and this seems like the norm based on what I've read about a couple of statistics packages. These values for $\epsilon$ are usually converted back into units of the original data, so the may appear differently in the R output.
Granger et al. (2010) points out that some authors have reported the value of the test statistic for a number of $\epsilon$ values and that this is not particularly helpful. My reading of the literature tells me that it is, indeed, the norm to report results using a number of values for $\epsilon$.
To test for non-linearity in practice, we must first build a linear model for a time-series $y_{t}$. For instance, we could follow the Box-Jenkins methodology and construct an ARIMA model. The BDS test is then applied to the residuals of this model. If the null hypothesis is rejected, one can conclude that there is non-linearity in the data, but the form of non-linearity cannot be determined by the test. Therefore, if non-linearity is present, we can't say whether it is chaos or a nonlinear stochastic process.
A second method of testing for non-linearity using the BDS test involves taking differences of the logs of the time-series and then applying the BDS test to this transformed series. For an example of this method, see Chu, P.K.K. (2001).
I've decided that this post is long enough, so I will pick up the application in R in the next post. Before wrapping up, though, a couple of final comments about the BDS test.
Granger et al. (2010) lend the following advice: at least in small samples up to 200 observations, the BDS test may not be the ideal choice. Interestingly, these course notes by Eric Zivot say that 500 observations is usually thought to be the minimal sample size for the BDS test to have reliable performance. Finally, one advantage of the BDS test is that it does not require any distributional assumption on the data being tested.
(Unlinked) Reference:
Brock, W., & Dechert, W. D. J. Scheinkman, 1987,“A Test for Independence Based on the Correlation Dimension.”. Economics Working Paper.
In this post I'd like to give a brief overview of the BDS test and at the end I'll provide an application of carrying out the BDS test in R. To introduce the BDS test, let me begin by saying that the method of the BDS test has its roots in the concept of spatial correlation from chaos theory.
Let's assume that $\{y_{t}\}$ is a strictly stationary process. Now consider the following vector called the m-history of the process: \[{\bf y}_{t}^{m} = (y_{t}, y_{t-1}, \ldots, y_{t+m-1})\] The construction of the m-history means that the univariate series $y_{t}$ has been embedded into a higher dimensional space. The purpose of this operation is to be able to examine the spatial correlation of the time-series. Notice that the vector, ${\bf y}_{t}^{m}$, consists of overlapping series. This becomes clear if we write out the following equations: \begin{aligned} y_{1}^{m} &= (y_{1}, y_{2}, \ldots, y_{m})\\ y_{2}^{m} &= (y_{2}, y_{3}, \ldots, y_{m+1})\\ & ~\vdots\\ y_{T-m+1}^{m} &= (y_{T-m}, y_{T-m+1}, \ldots, y_{T}) \end{aligned} To demonstrate the concept of embedding (which is due to Dutch mathematician Floris Takens), here is an example from Cromwell et al. (1994). Suppose we have a time-series, $y_{t}$, where $t = 1, 2, \ldots, 5$ and $m = 3$. Embedding this series into higher dimensions results in the following sequence of overlapping series: \begin{aligned} y_{1}^{3} &= (y_{1}, y_{2}, y_{3})\\ y_{2}^{3} &= (y_{2}, y_{3}, y_{4})\\ y_{3}^{3} &= (y_{3}, y_{4}, y_{5}) \end{aligned} In this example, the embedding operation creates three new three dimensional vectors: $y_{1}^{3}, y_{2}^{3}, y_{3}^{3}$. A requirement is that the vectors have the same length. Once the series has been embedded, we need an instrument to allows us to measure (a metric).
The correlation integral, introduced by Grassberger and Procaccia (1983), is a measure of spatial correlation. It measures the distance between the points in the m-dimensional space constructed by the embedding procedure; in the example above, it measures the distance between points in three dimensions. The correlation integral can also be interpreted as a measure of the frequency with which temporal patterns are repeated in the data.
The correlation integral at embedding dimension m is defined as: \[ C_{m, T} (\epsilon) = T_{m} (T_{m} - 1)^{-1} \sum_{t \neq s} I[y_{t}^{m}, y_{s}^{m}, \epsilon ] \] where $T_{m} = T - m + 1$, $t$ and $s$ both range from $1$ to $T - m + 1$ in the summation, m is the embedding dimension, and $\epsilon$ is the threshold value. $I(\cdot)$ is called the indicator function and is defined as: \begin{aligned} I[y_{t}^{m}, y_{s}^{m}, \epsilon ] &= 1 ~~~~~\text{if}~ || y_{t}^{m} - y_{s}^{m} || \leq \epsilon \\ &= 0 ~~~~~\text{otherwise} \end{aligned} The term, $|| y_{t}^{m} - y_{s}^{m} ||$, is a metric for the distance between two pairs. An example using the aforementioned hypothetical data set would be that the distance between $y_{1}^{3} = [y_{1}, y_{2}, y_{3}]$ and $y_{3}^{3} = [y_{3}, y_{4}, y_{5}]$ is given by: \[y_{1}^{3} - y_{3}^{3} = [ y_{1} - y_{3}, y_{2} - y_{4}, y_{3} - y_{5} ] \] The next step of the test involves choosing the maximum element from this resulting vector and then comparing it to the threshold value, denoted, $\epsilon$. If the maximum value is greater than $\epsilon$ then the pair is counted. Otherwise, the pair is not counted. The counting continues as the values for $t$ and $s$ are varied. Then, the number of counted pairs is divided by the total number of possible pairs. This fraction of pairs that are within a distance of $\epsilon$ of each other is the computed value for the correlation integral. The extreme cases are that all pairs are chosen, $C_{m, T} (\epsilon) = 1$, and that no pairs are chosen, $C_{m, T} (\epsilon) = 0$.
The BDS statistic is: \[ BDS = \sqrt{T} \frac{C_{m,T}(\epsilon) - C_{1,T}^{m}(\epsilon)}{\sigma_{m}(\epsilon)}\] where the variance, $\sigma_{m}^{2}(\epsilon)$, has the following complicated definition: \[ \sigma_{m}^{2}(\epsilon) = 4\{K(\epsilon)^{m} + 2 [\sum_{i=1}^{m-1} K(\epsilon)^{m-i} C(\epsilon)^{2i} + (m-1)^{2} C(\epsilon)^{2m} - m^{2} K (\epsilon) C(\epsilon)^{2m-2}]\} \] The indicator functions $C(\epsilon) = E\{I[y_{i}, y_{j}, \epsilon ]\}$ and $K(\epsilon) = E\{I[y_{i}, y_{j}, \epsilon ] I[y_{j}, y_{k}, \epsilon ]\}$ are defined in a similar way. For a complete definition, see Brock et al. (1996).
The null hypothesis of the BDS test is $H_{0}: y_{t}$ is i.i.d. The test is a two tailed test, so we should reject the null hypothesis of i.i.d at the 5% significance level when $|BDS| > 1.96$.
Notice that the statistic is a function of two parameters; the embedding dimension, m, and the threshold value, $\epsilon$. Values for these two parameters must be chosen, so this raises the question of how to choose values for m and $\epsilon$.
Various rules of thumb have been used to choose a value or values for $\epsilon$ and a short list of references can be found be referring to the concluding section of Brock (1995). The default values using the bds.test() function in R are $\epsilon$ = (0.5, 1.0, 1.5, 2.0) and this seems like the norm based on what I've read about a couple of statistics packages. These values for $\epsilon$ are usually converted back into units of the original data, so the may appear differently in the R output.
Granger et al. (2010) points out that some authors have reported the value of the test statistic for a number of $\epsilon$ values and that this is not particularly helpful. My reading of the literature tells me that it is, indeed, the norm to report results using a number of values for $\epsilon$.
To test for non-linearity in practice, we must first build a linear model for a time-series $y_{t}$. For instance, we could follow the Box-Jenkins methodology and construct an ARIMA model. The BDS test is then applied to the residuals of this model. If the null hypothesis is rejected, one can conclude that there is non-linearity in the data, but the form of non-linearity cannot be determined by the test. Therefore, if non-linearity is present, we can't say whether it is chaos or a nonlinear stochastic process.
A second method of testing for non-linearity using the BDS test involves taking differences of the logs of the time-series and then applying the BDS test to this transformed series. For an example of this method, see Chu, P.K.K. (2001).
I've decided that this post is long enough, so I will pick up the application in R in the next post. Before wrapping up, though, a couple of final comments about the BDS test.
Granger et al. (2010) lend the following advice: at least in small samples up to 200 observations, the BDS test may not be the ideal choice. Interestingly, these course notes by Eric Zivot say that 500 observations is usually thought to be the minimal sample size for the BDS test to have reliable performance. Finally, one advantage of the BDS test is that it does not require any distributional assumption on the data being tested.
(Unlinked) Reference:
Brock, W., & Dechert, W. D. J. Scheinkman, 1987,“A Test for Independence Based on the Correlation Dimension.”. Economics Working Paper.
Thursday, 2 May 2013
Exponentially Weighted Moving Average models
Open up any finance textbook and you won't have to look far before you come across exponentially weighted average (EWMA) models. The EWMA is typically used to measure historical volatility, but it can, of course, be used in other contexts. For example, see Muth (1960) for an early discussion of the statistical properties of EWMA models and its use in macroeconomics, in particular, in relation to the adaptive expectations mechanism associated with Milton Friedman.
Based on the amount of online video tutorials and other online content providing explanations of the EWMA, I can only infer that the EWMA is a popular tool for those working in the finance industry. And why not, eh? It's a simple model that can be understood and interpreted by mangers and the like.
Indeed, for some time-series, the EWMA can provide an adequate representation of the data generating mechanism and produce optimal forecasts, however, the EWMA can also produce sub-optimal forecasts when (mis)applied to other series of data. For this reason, using the EWMA simply out of habit may be a bad habit worth kicking.
Before using the EWMA for forecasting, one ought to grant some time to considering whether or not the EWMA is appropriate for the data in question. That is, of course, under the assumption that time is not of the essence and that the practitioners' ultimate goal is to produce optimal forecasts. In other words, let's assume that the practitioner is not under a severe time constraint and faced with the task of updating or developing hundreds of forecasting models for a company.
The claim thus far is that strict use of the EWMA is not the highest standard of forecasting. A justification of this claim now follows.
Consider the following fact:
Now, recall the fact that properly constructed ARIMA models produce optimal univariate forecasts - meaning that the mean-squared forecast error is smaller than for any other univariate forecast (produced by univariate, linear, fixed coefficient models).
Suppose now that rather than using the EWMA by habit, we instead use the Box-Jenkins method and let the data talk to us. In this situation, there are two possible outcomes; the appropriate model may, or may not, be an ARIMA(0,1,1), or, equivalently, an EWMA.
If it turns out that the appropriate model is an ARIMA(0,1,1) then our forecasts will be optimal. On the other hand, if the appropriate model is not an ARIMA(0,1,1) then our model has not been properly constructed and thus it must be admitted that our forecasts could be improved by specifying a different model from the suite of ARIMA models.
Based on the amount of online video tutorials and other online content providing explanations of the EWMA, I can only infer that the EWMA is a popular tool for those working in the finance industry. And why not, eh? It's a simple model that can be understood and interpreted by mangers and the like.
Indeed, for some time-series, the EWMA can provide an adequate representation of the data generating mechanism and produce optimal forecasts, however, the EWMA can also produce sub-optimal forecasts when (mis)applied to other series of data. For this reason, using the EWMA simply out of habit may be a bad habit worth kicking.
Before using the EWMA for forecasting, one ought to grant some time to considering whether or not the EWMA is appropriate for the data in question. That is, of course, under the assumption that time is not of the essence and that the practitioners' ultimate goal is to produce optimal forecasts. In other words, let's assume that the practitioner is not under a severe time constraint and faced with the task of updating or developing hundreds of forecasting models for a company.
The claim thus far is that strict use of the EWMA is not the highest standard of forecasting. A justification of this claim now follows.
Consider the following fact:
- The exponentially weighted moving average is algebraically equivalent to an ARIMA(0,1,1) model.
To put it another way, the EWMA is a particular model within the class of ARIMA models. In fact, there are various types of EWMA models and these happen to be included in the class of ARIMA(0,d,q) models - see Cogger (1974).
Now, recall the fact that properly constructed ARIMA models produce optimal univariate forecasts - meaning that the mean-squared forecast error is smaller than for any other univariate forecast (produced by univariate, linear, fixed coefficient models).
Suppose now that rather than using the EWMA by habit, we instead use the Box-Jenkins method and let the data talk to us. In this situation, there are two possible outcomes; the appropriate model may, or may not, be an ARIMA(0,1,1), or, equivalently, an EWMA.
If it turns out that the appropriate model is an ARIMA(0,1,1) then our forecasts will be optimal. On the other hand, if the appropriate model is not an ARIMA(0,1,1) then our model has not been properly constructed and thus it must be admitted that our forecasts could be improved by specifying a different model from the suite of ARIMA models.
R Time Series Issue: Seasonal ARIMA Models
In my previous post, I mentioned that there are seasonal non-seasonal (S-NS) multiplicative models and seasonal non-seasonal additive models.
About these two types of models, Pankratz (1983) tells us that:
To fit a S-NS model in R, we can specify the function using code in the following form:
I'm not the first to raise this question about S-NS models in R. The question was asked on the R mailing list here, but no answer was given. A follow up question was asked and in the email it is revealed that the arima() function assumes the models to be multiplicative.
The question now is: how to model S-NS additive models in R?
For the moment, one way around the problem of estimating S-NS additive models would be to use the lm() function. A proper solution would be to write code to slightly revamp arima(). I'm not sure if anyone is working on this problem, so I might try and write my own function that will do the trick.
About these two types of models, Pankratz (1983) tells us that:
Multiplicative models are an excellent starting place for modeling S-NS realizations, and they often give good results. Most computer programs for estimating UBJ-ARIMA models provide at least an option for multiplicative seasonality; some simply assume that seasonal models are multiplicative.It is then said that:
Occasionally, the assumption that the seasonal and nonseasonal elements are multiplicative is inappropriate.This led me to question whether or not S-NS models in R using the arima() function are multiplicative by assumption or if there is an option to choose between multiplicative and additive models.
To fit a S-NS model in R, we can specify the function using code in the following form:
arima(x, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = s)
Other options can be specified, but there doesn't appear to be any reference to say whether or not the models are multiplicative by assumption.I'm not the first to raise this question about S-NS models in R. The question was asked on the R mailing list here, but no answer was given. A follow up question was asked and in the email it is revealed that the arima() function assumes the models to be multiplicative.
The question now is: how to model S-NS additive models in R?
For the moment, one way around the problem of estimating S-NS additive models would be to use the lm() function. A proper solution would be to write code to slightly revamp arima(). I'm not sure if anyone is working on this problem, so I might try and write my own function that will do the trick.
Subscribe to:
Posts (Atom)