# Estimating accurately the mean of an autocorrelated bounded integer time series

I have a bounded integer time series $X_{1:\infty}$ ($1\leq X_k\leq M$), and I want to estimate the mean $$s = \lim_{L\to\infty} \frac{1}{L}\sum_{k=1}^L X_k.$$ I'm assuming it exists and that $X_k$ converges in distribution to some distribution with mean $\mu$ and variance $\sigma^2$, and $s$ converges to $\mu$.

I am estimating it for $X_{1:L}$ as $s_X = L^{-1}\sum X_{1:L}$, and the accuracy of $s_X$ I can estimate using $$\mathbb{V}[s_X] = \frac{\sigma^2}{L} + \frac{2}{\sigma^2L} \sum_{\tau\geq1} \mathrm{Acf}(\tau)(1-\tau/L),$$ by replacing $\sigma$ and the autocorrelation function $\mathrm{Acf}(\tau)$ with the sample standard deviation and the sample a.c.f..

The asymptotic distribution of $X_k$ seems to be quite unusual, with 75% of values being $\leq 4$, and approximately power-law-like distributed up to the maximum $M\sim10^4$. The mean and s.d. are about 12 and 150. The autocorrelation function seems to decay roughly exponentially as $\mathrm{Acf}(\tau)\sim 2^{-\tau}$ until trailing off at $\mathrm{Acf}(\tau)\approx 10^{-3}$.

So with these unusual properties of $X_k$, how well can I expect $s_X$ to estimate the limit, and how can I check how good it is?

More importantly, are there some better estimators that I can try?

• You may be able to derive an anylitical solution with the conditional distribution of $X_k$. I.e. $f (X_k | X_1,...,X_{k-1})$. Can you tell us what that conditional distribution is? Other than that you can always simulate the process and estimate the mean many times which will give you an emperical approximation for the estimators variance. The estimators bias could also be approximating with bootstrapping techniques. All that said, it is still probably better to work with an analytical solution using the conditions where possible. – Zachary Blumenfeld Feb 7 '16 at 0:19
• @ZacharyBlumenfeld I am unable to derive the distribution in closed form. Simulating it is what I've been doing so far: my question is what's the best way to extract an accurate estimate from the process sample I have. Is it possible to do better? – Kirill Feb 7 '16 at 15:31
• The answer is yes. You can use Feasible Generalized Least Squares (FGLS) to produce a more efficient estimate. This is a special enough case though that an answer can be made really simple. Will try to get to it later. – Zachary Blumenfeld Feb 8 '16 at 1:10
• You seem to be estimating the variance of the mean without smoothing. This will result in an inconsistent estimator of the variance. In such a case you would need to use "fixed-b" asymptotics. Two questions (I happen to be aware of because I contributed) dealing with the issue are stats.stackexchange.com/questions/153444/… and stats.stackexchange.com/questions/60942/… – Christoph Hanck Feb 9 '16 at 7:48

## The Problem

To re-state the problem: You have a time series $\{X_t\}_{t=1}^\infty$ and are interested in estimating the unconditional mean, $E[X_t]$.

You do not know the conditional distribution of $X_t$ (i.e. $f(X_t|X_1,...,X_{t-1})$ is unknown). However, you do know that the time series is stationary and ergodic (obviously) and that it is homoscedastic; $var(X_t)=\sigma^2\;\;\forall t \in \mathbb{N}$.

Also, for any lag $k \in \mathbb{N}$ you know the auto-correlation coefficient $\rho_k = cor(X_t,X_{t-k})$ (at least you can approximate it well in a large enough sample). With this information, you can write the covariance matrix for the finite sample $X_1,...,X_T$; $$\mathbf{cov}(X_1,...,X_T)=\Sigma = \sigma^2 R$$ where $$R=\begin{bmatrix} 1 & \rho_{1} & \cdots & \rho_{T-1} \\ \rho_{1} & 1 & \cdots & \rho_{T-2} \\ \vdots & \vdots & \ddots & \vdots \\ \rho_{T-1} & \rho_{T-2} & \cdots & 1 \end{bmatrix}$$

## Deriving The Best Linear Unbiased Estimator (BLUE)

(See the Wikipedia entry on GLS and the references cited therein for more details on what follows)

Define the vectors $\mathbf{X} = (X_1,X_2,...,X_T)$ and $\boldsymbol{\mu}=(\mu,...,\mu) \in \mathbb{R}^T$. You can consistently estimate $E[X_t]$, by minimizing the following quadratic $$\hat\mu_{W} = \min_\mu\bigg\{(\mathbf{X}-\boldsymbol{\mu})'W(\mathbf{X}-\boldsymbol{\mu}) \bigg\}$$ Where $W$ is a $T\times T$ positive-definite symmetric matrix. The closed form solution is $$\hat\mu_{W} = (\mathbf{1}'W\mathbf{1})^{-1}\mathbf{1}'W\mathbf{X}$$

where $\mathbf{1}$ is a $T\times 1$ vector of ones. The common sample average $\frac{1}{T}\sum_{t=1}^T X_t$ is a special case of the above where we choose $W \propto I$, where $I$ is the identity matrix*.

*As a side note, $\hat\mu_{W}$ is indifferent to proportionality constants in $W$. i.e. $\forall\; a \in \mathbb{R}_{++}$, $\hat\mu_{W}=\hat\mu_{aW}$.

Any choice of $W$, so long as it is positive-definite and symmetric, will result in a consistent estimator of $E[X_t]$. In other-words, it can be shown that as $T\rightarrow \infty$; $\hat\mu_{W}\stackrel{p}{\rightarrow} E[X_t]$, $\forall\; W \in \mathcal{S}$,where $\mathcal{S} \subset \mathbb{R}^{T\times T}$ is the subset of all positive definite symmetric matrices.

However, some choices of $W$ yield more efficient (lower variance) estimates than others. It can be shown that the most efficient choice of $W$, resulting in the Best-Linear Unbiased Estimate (BLUE) of $E[X_t]$, is $W \propto \Sigma^{-1}$. In your case, you can use;

$$\boxed{ \hat\mu_{R^{-1}} = (\mathbf{1}'R^{-1}\mathbf{1})^{-1}\mathbf{1}'R^{-1}\mathbf{X} }$$

## Variance and Other Asymptotic Properties

$\hat\mu_{R^{-1}}$ is unbiased, consistent, efficient, and asymptotically normal: $$(\hat\mu_{R^{-1}} - E[X_t])\ \xrightarrow{d}\ \mathcal{N}\!\left(0,\,\sigma^2(\mathbf{1}'R^{-1}\mathbf{1})^{-1}\right)$$ So $$var(\hat\mu_{R^{-1}}) \approx \sigma^2(\mathbf{1}'R^{-1}\mathbf{1})^{-1}$$

Also, even though $\hat\mu_{R^{-1}}$ has a lower variance than $\hat\mu_{I}=\frac{1}{T}\sum_{t=1}^T X_t$, the difference may be negligibly small. If that is the case you may prefer to use $\frac{1}{T}\sum_{t=1}^T X_t$ since calculating $\hat\mu_{R^{-1}}$ is much more computationally expensive. You can calculate the difference in variance analytically to get an idea.
1. simulate $(X_1,...X_T)$ and collect $\hat\mu_{R^{-1}}$ and $\hat\mu_{I}$ multiple times. If you are given the data $X_1,...X_T$ and cannot directly simulate it, you can use a block-bootstrap.
2. Plot the histogram of $\hat\mu_{R^{-1}}$ and $\hat\mu_{I}$. Do they look approximately normally distributed? Are the sample variance of the estimates close to the asymptotic approximations?
If the asymptotics don't seem to approximate well, you could either increase $T$ or use the empirical results of the simulation in place of the asymptotic results.