# OLS estimation of intercept in AR($p$) in R

I investigate the performance of the OLS estimator of an AR($$3$$) model given by $$X_t=\mu+\phi_1X_{t-1}+\phi_2X_{t-2}+\phi_3X_{t-3}+\varepsilon_t$$ for $$t\in\mathbb Z$$ using the following code:

set.seed(123)
est<-matrix(NA,nrow=10^3,ncol=5)
for (j in 1:10^3)
{
x<-5+arima.sim(list(ar=c(.4,.3,.2)),10^3)
m<-ar.ols(x,aic=FALSE,order.max=3,demean=FALSE,intercept=TRUE)
est[j,1]<-m$$x.intercept est[j,2:4]<-m$$ar
est[j,5]<-m$var.pred } colMeans(est)  These are the results that I get: [1] 0.5377802 0.3983299 0.2967576 0.1971440 0.9954128  All of the parameters are estimated reasonably well except $$\mu$$. It seems that $$0.1\mu$$ is estimated instead of $$\mu$$ since $$\mu=5$$. In the description of ar.ols, it is stated that: Some care is needed if intercept is true and demean is false. Only use this if the series are roughly centred on zero. Otherwise the computations may be inaccurate or fail entirely. $$\{X_t\}_{t\in\mathbb Z}$$ is not centred around zero so this is probably where the problem is. However, I do not understand why this is the case. The OLS estimator is consistent in this case and I think that $$\mu$$ should also be estimated reasonably well. Could someone explain where the problem is and how I should estimate $$\mu$$ correctly? Any help is much appreciated! ## 1 Answer What you simulate is (the arima.sim part) $$Y_t=\phi_1Y_{t-1}+\phi_2Y_{t-2}+\phi_3Y_{t-3}+\varepsilon_t$$ and, for m <- 5, $$X_t=\mu+Y_t$$, or $$X_t-m=\phi_1(X_{t-1}-m)+\phi_2(X_{t-2}-m)+\phi_3(X_{t-3}-m)+\varepsilon_t$$ or $$X_t=m\left(1-\sum_{j=1}^3\phi_j\right)+\phi_1X_{t-1}+\phi_2X_{t-2}+\phi_3X_{t-3}+\varepsilon_t$$ where $$m\left(1-\sum_{j=1}^3\phi_j\right)=5(1-0.9)=0.5$$, which is actually estimated pretty well, with an average of 0.5377802. This is therefore not the same as simulating $$X_t=m+\phi_1X_{t-1}+\phi_2X_{t-2}+\phi_3X_{t-3}+\varepsilon_t$$ The issue is related to the fact that the expected value of an AR(p) process $$Y_t=c+\phi_1Y_{t-1}+\phi_2Y_{t-2}+\ldots+\phi_pY_{t-p}+\varepsilon_t$$ is given by (see e.g. here) $$\mu=\frac{c}{1-\sum_{j=1}^p\phi_j}$$ When taking m even larger, computations may fail entirely, as mentioned by ?ar.ols. Taking an AR(1) for simplicity, we get > n <- 10^2 > x <- 500 + arima.sim(list(ar=.4), n) > m <- ar.ols(x, aic=FALSE, order.max=1, demean=FALSE, intercept=TRUE) Error in AA %*% t(X) : requires numeric/complex matrix/vector arguments In addition: Warning message: In ar.ols(x, aic = FALSE, order.max = 1, demean = FALSE, intercept = TRUE) : model order: 1 singularities in the computation of the projection matrix results are only valid up to model order 0  The "culprit" in ar.ols seems to be in the lines where the regressor matrix is generated: # example data n <- 10^2 c.int <- 500 x <- c.int + arima.sim(list(ar=.4), n) # taken & simplified from the source code: m <- 1 # the single lag y <- embed(x, m + 1L) # generates lagged y X <- cbind(rep.int(1, nrow(y)), y[, 2:ncol(y)]) # adds the constant for the intercept XX <- t(X) %*% X rank <- qr(XX)$rank
if (rank != nrow(XX)) {
warning(paste("model order: ", m, "singularities in the computation of the projection matrix",
"results are only valid up to model order", m - 1L), domain = NA)
}


Output:

> XX
[,1]        [,2]
[1,]    99.00    49491.75
[2,] 49491.75 24741893.87

> rank
[1] 1


Effectively, the sum of squares of the regressors (the 2,2-element of XX) seems to be too large relative to n (the 1,1-element plus 1) to ensure, according to ar.ols, sufficiently high numerical precision of the $$(X'X)^{-1}$$ matrix in the OLS estimate.

Interestingly, other routines report a rank of 2 for this matrix, corresponding to what we know to apply in this example:

> library(Matrix)

> rankMatrix(XX)
[1] 2

• Thank you very much for such a thorough answer! Your answer also explains why some care is needed if intercept is true and demean is false in the code, right? The thing is that the intercept is not equal to the mean unless the mean is zero. Is that correct? – Cm7F7Bb Mar 7 '19 at 10:42
• Yes, if $c\approx0$, then so is the mean of the process and the distinction as spelled out in the answer does not matter. You may also find this thread to be useful: stats.stackexchange.com/questions/176208/… – Christoph Hanck Mar 7 '19 at 10:44
• I still don't understand why we can only set intercept to true and demean to false if the series are roughly centred on zero. If we do not demean the data, we have have an AR($p$) model with an intercept and the OLS estimator is consistent in this case. The simulation that I have shows that this works and I don't observe any problems. Why do they say that the computations may be inaccurate or fail entirely? Would you have an explanation? – Cm7F7Bb Mar 7 '19 at 11:26
• That is not so clear indeed. One could dig in how ar.ols precisely proceeds. Indeed, n <- 10^4 m <- 500 x <- m + arima.sim(list(ar=.4), n) summary(lm(x[2:n]~x[1:(n-1)])) ar.ols(x, aic=FALSE, order.max=1, demean=FALSE, intercept=TRUE) reveals that using lm works correctly. – Christoph Hanck Mar 7 '19 at 11:47
• I made another edit which gives at least some indications. – Christoph Hanck Mar 7 '19 at 12:36