3
$\begingroup$

I want to estimate a vector-valued model

$$\mathbf{y}_t = a\mathbf{y}_{t-1}+b\mathbf{y}_{t-2}+\cdots$$

Here, each $\mathbf{y}_t\in\mathbb{R}^n$ and the coefficients $a,b,\dotsc$ are real numbers (unlike in VAR, where they are matrices). I wonder what the best way to do this in R is. Applying the function ar doesn't seem to work for me ([i] is this only for series that lie in $\mathbb{R}^1$?), and when I apply standard linear regression, I'm not sure how I should test for the order of the model, etc. [ii] Should I go for the $R^2$ value and significance of the coefficients?

Also, I have my data in a matrix $Y$ where each column corresponds to one time point. Currently, I regress as follows

Y = c(Y[,10], Y[,9],..., Y[,2])
Y.l1 = c(Y[,9],..., Y[,1])
fit = lm(Y ~ 0 + Y.l1)

[iii] Is there a more elegant way?

$\endgroup$
  • $\begingroup$ You have at least three different questions in the post (How to technically implement the model? Does ar work with univariate data only? How to select model order?), but they seem to get lost in the flow of the text. You may consider highlighting them by restructuring your post. $\endgroup$ – Richard Hardy Nov 27 '15 at 10:41
3
$\begingroup$

First of all let us state that you are trying to fit the model which is not that common. Since $y_t\in \mathbb{R}^n$ you can always estimate $n$ independent AR models. So your model implies that all these $n$ AR models must have the same coefficients. That is a powerful restriction and you need to be sure, that the process you are trying to model satisfies this restriction.

Having said that, if you are sure that the restriction is feasible, the model is a linear model, so you can use OLS to fit it. Naturally R function ar is not suitable for this model, because as I said such model is not that common, so nobody thought that it is worth implementing.

The model order can be chosen in the same way VAR model order is chosen, i.e. using information criteria. If you fit the model via OLS, you can use AIC and BIC with no problems. Just remember that you need to compare AIC and BIC for models which are fitted on the same sample. For VAR models this means dropping the first $k$ observations, where $k$ is the number of highest lag number you are trying to estimate.

For elegant way of fitting (for my understanding of elegant), I recommend using dplyr and tidyr packages to manipulate your data and then fit using lm.

Suppose your data is a matrix, where data for time periods are in columns:

set.seed(100)
Y<-matrix(rnorm(1000),nrow = 10)
colnames(Y) <- 1:100

> Y[1:5, 1:2]

               1           2
[1,] -0.50219235  0.08988614
[2,]  0.13153117  0.09627446
[3,] -0.07891709 -0.20163395
[4,]  0.88678481  0.73984050
[5,]  0.11697127  0.12337950

Your data is in a wide format and I suggest it to convert to long format:

library(dplyr)
library(tidyr)

dY <- Y %>% 
data.frame(check.names = FALSE) %>% 
mutate(id = 1:10)  %>% 
gather("time",Y,-id, convert = TRUE)

> dY[1:10,]
   id time       Y
1   1    1 -0.50219235
2   2    1  0.13153117
3   3    1 -0.07891709
4   4    1  0.88678481
5   5    1  0.11697127
6   6    1  0.31863009
7   7    1 -0.58179068
8   8    1  0.71453271
9   9    1 -0.82525943
10 10    1 -0.35986213

In this format you have two additional columns which record the time period and to which element of the vector the value corresponds to.

Having such format it is easy to add lagged values of $y$:

dY <- dY %>% group_by(id) %>% mutate(Y.l1 = lag(Y,1), Y.l2 = lag(Y, 2), Y.l3 = lag(Y,3))

To fit the model use lm:

lm(Y~Y.l1,data= dY)

Call:
lm(formula = Y ~ Y.l1, data = dY)

Coefficients:
(Intercept)         Y.l1  
    0.01749     -0.01624  

To choose the order, you need to fit the model for all orders:

ar1 <- lm(Y~Y.l1, data = dY %>% filter(time > 3))
ar2 <- lm(Y~Y.l1+Y.l2, data = dY %>% filter(time > 3))
ar3 <- lm(Y~Y.l1+Y.l2+Y.l3, data = dY %>% filter(time > 3))

Put them all in the list for easier manipulation:

mods <- list(ar1, ar2, ar3)

Check that the sample size is the same for all models, before calculating AIC/BIC:

> mods %>% sapply(function(x)nrow(x$model))
[1] 970 970 970

Calculate AIC:

> mods %>% sapply(AIC)
[1] 2832.232 2832.308 2832.659

Calculate BIC:

> mods %>% sapply(BIC)
[1] 2846.864 2851.817 2857.045
$\endgroup$
  • $\begingroup$ Thanks a lot mpiktas. I have to digest this now but I think your answer is very helpful (writing from a newly generated account) $\endgroup$ – mathse Nov 27 '15 at 13:02
  • $\begingroup$ If you have created two accounts you can ask to merge them: stats.stackexchange.com/help/merging-accounts $\endgroup$ – mpiktas Nov 27 '15 at 14:21
  • $\begingroup$ Hi @mpiktas, I messaged them. But they seem to ignore it, although they should easily have seen my IP. So, I'm sorry that I can't accept your answer. $\endgroup$ – mathse Dec 1 '15 at 13:18
  • $\begingroup$ No worries. Maybe moderators are busy at the moment. $\endgroup$ – mpiktas Dec 1 '15 at 13:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.