Re: AR model

[Juan Pablo Carbajal]

> > I have a doubt on polynomials. If I have an AR model like
> > Yt=b0+b1*Yt-1+b2*Yt-2+et and I want to write the corresponding
> > polynomial, would it be [1 -b1 -b2]?
> You seems to use negative ones b0+b1*Yt^-1+b2*Yt^-2
> so they are not polynomial

The polynomial representation is not the delayed representation, in estefania's
x^{-2} means x(t-2*dt), that is the value of the signals two timesteps
in the past. This is a operator representation. the corresponding
polynomial would be
x = a x^{-2} --> x(t) = a x(t-2*dt), so the coefficients are the same.

@estefania:
Yt=b0+b1*Yt-1+b2*Yt-2+et
is then
y(t) = b0 + b1 * y(t-1) + b2 * y(t-2) + noise(t)
and the polynomial is
py = [b2 b1 b0]

if you want to solve the coefficients that minimize the squared error
(i.e. is the same as assuming noise is i.i.d. and sampled from a
Gaussian distribution) you can do
(for one data set)
# assumes x is a column vector containing the signal, with length N
X = [x(1:end-2) x(2:end-1) ones(N-2,1)]; # this is the Vandermonde
matrix in operator space
py = X \ x(3:N)

Example:

# create synthetic data
N = 10;
x_noiseless = x = zeros (N, 1);
x_noiseless(1) = 0.0;
x_noiseless(2) = 0.0;
px_true = [0.2 -0.5 0.1].';
pdeg = length(px_true) - 1; # polynomial's degree
for i = 3:N
    x_noiseless(i) = px_true(3) + x_noiseless(i-1) * px_true(2) +
x_noiseless(i-2) * px_true(1);
    # alternative, and more generic
    # x_noiseless(i) = [x((i-pdeg):(i-1)).' 1] * px_true;
endfor
x = x_noiseless + 0.01 * randn(N, 1); # add Gaussian noise

# Solve with least-squares (check other functions in optim package,
and other optimizations functions)
X = [x(1:N-2) x(2:N-1) ones(N-2,1)];
px = X \ x(3:N);

printf ('px_true and px_ols\n');
disp([px_true px])

# Compare
x_ = zeros(N, 1);
for i =3:N
    x_(i) = [x((i-pdeg):(i-1)).' 1] * px;
endfor
plot(x,'o;data;', x_, 'x;estimation: mean value;')

# end example

You can vectorize the for loop using the function movslice, or if you
implement your Ar model as a function you can use movfun.
Check their help

Regards,