How many samples are needed for guaranteed maximum estimation error of the mean when using sample mean estimator?

Question

I am estimating the mean of a Linear Quadratic Regulation process(LQR) using the sample mean estimator. It seems to me that my estimation error is large with respect to the number of samples that I use. Therefore, I am wondering if there exists some theory about the size of the estimation error with respect to the number of samples.

Details:

Let $x_t\in\mathbb{R}$, $u_t\in\mathbb{R}$ and $v_t\in\mathbb{R}$ be random processes such that \begin{equation} \begin{split} x_{t+1} &= Ax_{t} + Bu_t + v_t\\ u_t & = -Kx_t \\ v_t & \sim \mathcal{N}(0, \sigma_v^2) \end{split} \end{equation} where $K$ is the LQR feedback gain.

Using this model I generate $N$ trajectories and calculate the sample mean of $x_t$. Let $x_{t,i}$ denote the $i$:th sample then the sample mean is given by \begin{equation} \bar{x}_t = \frac{1}{N} \sum_{i=1}^N x_{t,i}. \end{equation}

The analytical mean of $x_t$ can be shown to be $\hat{x}_t = (A-BK)\hat{x}_{t-1}$. Thus it can be calculated given an initial condition.

Computing $\hat{x}_t - \bar{x}_t$ for $N = 10^7$ I get an error sizing $10^{-4}$, which I believe is too large with respect to the number of samples that I use.

Here is my Matlab code:

clear all;

A = 1.05;
B = 1;
sigma2_v = 0.5^2;
Q = 1.0;
R = 1.0;
N = 0.0;
K = dlqr(A,B,Q,R,N);

nbr_trajectories = 10000000;
time_horizon = 6;
initial_state = 10.0;

% Generate data
trajectories = cell(nbr_trajectories, 1);

for i = 1 : nbr_trajectories
    trajectory = zeros(time_horizon, 1);
    state = initial_state;
    for j = 1 : time_horizon
        trajectory(j) = state;
        state = (A-B*K)*state + sqrt(sigma2_v) * randn();
    end
    trajectories{i} = trajectory;
end

% Estimate sample mean
sample_mean = zeros(time_horizon, 1);
for i = 1 : time_horizon

    % Estimate mean
    mean = 0.0;
    for j = 1 : nbr_trajectories
        mean = mean + trajectories{j}(i);
    end
    mean = mean / nbr_trajectories;
    sample_mean(i) = mean;
end

% Calculate analytic mean
analytic_mean = zeros(time_horizon, 1);
state_mean = initial_state;
for i = 1 : time_horizon

    analytic_mean(i) = state_mean;

    % Update next timestep;
    state_mean = (A-B*K) * state_mean;
end

% Compare sample mean and analytic mean
for i = 1 : time_horizon - 1
    strcat('i = ', num2str(i))
    sample_mean(i) - analytic_mean(i)
end