
From:  estefaniame 
Subject:  Re: AR model 
Date:  Sun, 5 Jul 2020 21:39:23 +0200 
> > I have a doubt on polynomials. If I have an AR model like
> > Yt=b0+b1*Yt1+b2*Yt2+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(t2*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(t2*dt), so the coefficients are the same.
@estefania:
Yt=b0+b1*Yt1+b2*Yt2+et
is then
y(t) = b0 + b1 * y(t1) + b2 * y(t2) + 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:end2) x(2:end1) ones(N2,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(i1) * px_true(2) +
x_noiseless(i2) * px_true(1);
# alternative, and more generic
# x_noiseless(i) = [x((ipdeg):(i1)).' 1] * px_true;
endfor
x = x_noiseless + 0.01 * randn(N, 1); # add Gaussian noise
# Solve with leastsquares (check other functions in optim package,
and other optimizations functions)
X = [x(1:N2) x(2:N1) ones(N2,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((ipdeg):(i1)).' 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,
[Prev in Thread]  Current Thread  [Next in Thread] 