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!
 A: 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

