Matlab: How do I implement the Expectation operator for time series
I'm not sure I understand the question. If I'm missing something, leave a comment and I'll try to edit/expand/clarify as needed.

First of all, "mean" and "expected value" are related concepts, but are not the same thing. You find the mean of a set of samples, or observations. If you sampled/observed properly, the mean will be close to the expected value, which is the theoretical, "true mean" of the underlying random process.

When you sample a random process, you end up with a set of random numbers. If you sample the same process again, you end up with a different set of numbers. So, you can think of each sample as a random variable (a black box that produces random numbers); the random process is then a sequence of random variables.

These random variables may be correlated; that means, intuitively, that if one of them has a certain value, you can use it to predict the values of some of the other samples. For example, you may find that if one sample is positive, its neighboring samples tend to be positive too.

In your question you say that you have one set of samples, y(t). In general, you don't have enough information to generate the correlation matrix. However, you could do it like this:

M = length(y);
C = zeros(M,M);
for i = 1:M
    for j = 1:M
        C(i,j) = mean(y(i)*y(j));  % Note that the mean is useless here,
                                   % since y(i)*y(j) is scalar
    end
end


To find a closer approximation to the actual correlation matrix of the underlying process, you need to gather many observations. Let's say that you have stored the observations in the rows of a matrix Y (every row is a different set of observations). Then you could do:

M = size(Y,2);
C = zeros(M,M);
for i = 1:M
    for j = 1:M
        C(i,j) = mean(Y(:,i).*Y(:,j));  
    end
end


To test this, let's try sampling a Gaussian process of zero mean and variance 1. Every time I look at the process, I'll gather three samples, and I'll look at the process 1000 times:

Y = randn(1000,3);


If I run the code above with matrix Y, I obtain this correlation matrix:

0.9598    0.0590    0.0238
0.0590    0.9675   -0.0286
0.0238   -0.0286    1.0055


which is expected (the "true" correlation matrix is the identity).

Regarding calculating a time shift in Matlab, it's helpful to arrange vectors in a table as follows (I've assumed $t=1,2,\ldots$ and $y(t)=10t$):

t    | y(t)  |  y(t-1)  | y(t-2) 
     |       |          |
1    |  10   |  --      | --
2    |  20   |  10      | --
3    |  30   |  20      | 10
4    |  40   |  30      | 20
...


It's clear from the table that, to find yn=y(t-n), you cut n elements from the end of y, and shift the time vector:

tn = t(n:end);    % shifted time vector
yn = y(1:end-n);  % shifted signal


The following code provides an example, which doesn't assume integer time or integer delay:

t = 0:0.1:3;  % time vector
y = exp(t/5); % y(t)
n = 0.8;      % time delay
tn = t(t&gt;n);  % delayed time vector
yn = y(1:length(tn));  % delayed signal
plot(t,y,'r',tn,yn,'g');  % plot


This produces the following plot: