# ARMA Estimation Efficiency of the Mean

Usually in books related to ARMA Time Series, it is assumed that the series is 0 mean.

If not, the recommended standard procedure is to simply subtract the sample mean from the series and continue business as usual.

My question is whether there is a difference in estimation efficiency of ARMA parameters and the mean parameter when estimated altogether with MLE procedure vs first subtracting the sample mean and then estimating only ARMA parameters with MLE procedure.

Yes, there is, but in practice it may be meaninglessly small.

Consider an ARMA process with unknown mean as a subset of the more general problem of a non-diagonal, possibly heteroskedastic covariance matrix in a linear regression of the form:

$$y = \beta X + \epsilon, \,\,\epsilon \sim \text{Normal}(0,\Sigma)$$

In your case, $$X$$ consists of a single column of ones, i.e., the intercept, and $$\beta$$ is the mean. In ARMA modeling, $$\Sigma$$ has a lot of structure, and can be written in terms of a finite number (this is important for asymptotics) of parameters $$\theta$$ (for the AR part) and $$\phi$$ (for the MA part): $$\Sigma(\theta, \phi)$$.

To see how the parameters of an ARMA process define a covariance matrix, consider the autocorrelation function (ACF) commonly used to help identify the lag structure of $$\theta$$ and $$\phi$$. This function gives the correlation between $$\epsilon_t$$ and $$\epsilon_{t+h}$$ for all $$h$$; clearly the relationship between this function and the ARMA parameters is is sufficient to define the correlation matrix of $$\epsilon$$ in terms of $$\theta$$ and $$\phi$$, and is one small step from defining the covariance matrix of $$\epsilon$$ in terms of $$\theta$$ and $$\phi$$.

For example, the correlation matrix associated with an AR(1) process is defined by:

$$\rho_{ij} = \theta^{|i-j|}$$

and that of an MA(1) process by:

$$\rho_{ii} = 1,\,\rho_{i,i+1} = \rho_{i,i-1} = {-\phi \over 1+\phi^2}$$

with other entries equal to $$0$$.

The two-step process described in your question essentially estimates $$\beta$$ via OLS, then the parameters of $$\Sigma(\theta, \phi)$$ via, for example, MLE. However, the Gauss-Markov theorem tells us that the OLS estimate of $$\beta$$ is not efficient when $$\Sigma$$ is known, instead, $$\tilde{\beta} = (X^T\Sigma^{-1}X)^{-1}X^T\Sigma^{-1}y$$ is the BLUE (Best Linear Unbiased Estimator). This is known as the Generalized Least Squares (GLS) estimator. This is also the MLE under the assumption of Normality of $$\epsilon$$.

Since we have to estimate $$\Sigma$$, Gauss-Markov doesn't apply, and we have to resort to some form of Feasible Generalized Least Squares (FGLS.) However, it turns out that the FGLS estimator will have the same asymptotic distribution as the GLS estimator if:

$$\text{plim} {X^T(\hat{\Sigma}^{-1} - \Sigma^{-1})X \over N} = 0$$

and

$$\text{plim} {X^T(\hat{\Sigma}^{-1} - \Sigma^{-1})\epsilon \over \sqrt{N}} = 0$$

Given that there are a finite number of parameters in $$\Sigma$$, and $$X$$ is just a vector of ones, these are not particularly strenuous assumptions to make. These are sufficient conditions, and, clearly, if they hold, the joint estimator is superior asymptotically to the two-step estimator.

When doing ARMA modeling, you would just use the standard time series algorithms for estimating ARMA parameters, not try to estimate the covariance matrix directly.

In practical applications, though, we are often not in the happy asymptotic world, and it's not so clear in finite samples that there is anything to be gained. Here's an example comparing the two procedures on a zero-mean ARMA(2,2) process with 250 observations:

parms <- list(ar=c(0.5,-0.2), ma=c(0.4,0.2))

res_one_step <- matrix(0, 10000, 5)
res_two_step <- matrix(0, 10000, 5)
for (i in 1:10000) {
x <- arima.sim(parms, n=250)
res_one_step[i,] <- as.numeric(arima(x, order=c(2,0,2))$$coef) res_two_step[i,5] <- mean(x) res_two_step[i,1:4] <- as.numeric(arima(x-res_two_step[i,5], order=c(2,0,2), include.mean=FALSE)$$coef)
}

mse_one_step <- apply(res_one_step, 2, var) + (colMeans(res_one_step) - c(0.5,-0.2,0.4,0.2,0))^2
mse_two_step <- apply(res_two_step, 2, var) + (colMeans(res_two_step) - c(0.5,-0.2,0.4,0.2,0))^2


with results:

> sqrt(mse_one_step)
[1] 0.3933203 0.1974865 0.4008505 0.1942562 0.1434249
> sqrt(mse_two_step)
[1] 0.3936735 0.1977207 0.4011226 0.1942403 0.1434510
> sqrt(mse_two_step)/sqrt(mse_one_step)
[1] 1.0008980 1.0011859 1.0006788 0.9999183 1.0001815


Essentially no difference. The fact that we are estimating the parameters of the true process instead of going through an identification or model selection step is unlikely to have made a difference, given how similar the results are.

• So we are assuming a sufficiently large covariance matrix to accommodate AR parameters, I presume, as the correlation decays slowly. – Cagdas Ozgenc Dec 17 '19 at 17:07
• No, we are just assuming the covariance matrix is appropriately sized to the data. For example, if you have an AR(1) with $\theta = 0.9$, but only 10 observations, you can still construct a perfectly valid covariance matrix, even though $0.95^{10}$ isn't all that small. Any ARMA process with any sample size defines the covariance matrix. Of course, for asymptotics, your sample size goes to $\infty$, so the size of your covariance matrix goes to $\infty \times \infty$, but it is still a function of a finite (and fixed) number of parameters, so it's OK. – jbowman Dec 17 '19 at 17:32

Provided that it is a not a rule to subtract to sample mean from the series before fitting a model (de-meaning is possible but not a gold standard in all applications), then we could say the following.

Theoretically, you are perfectly right: all the parameters of MLE should be estimated altogether. So, in this case, also the sample mean, which is indeed the sample estimate of the unconditional mean. However, the real problem for the efficiency of the MLE comes when you make the WRONG assumption about the true value of the parameters that you have pre-estimated. In this case, if you have enough data, the estimation of the sample mean can be considered reliable as a proxy for the population mean. If this assumption is correct and you have correctly specified the model, then you should have no relevant loss of efficiency. Instead it is difficult to pre-estimate reliably other information (think about the shape of the distribution of innovations for example). That is a different pair of shoes, very “Error-prone”.