%
% serie-storiche:
%
N = 256;
clear i;
a = poly([0.3 0.7+0.4i 0.7-0.4i]);   a = a/a(1);
na = length(a);
y = 7.345*exp(0.0097*(0:N-1)) + filter([1],a,randn(1,N));
figure(3); clf, hold on, plot(0:N-1,y,'r-'); pause,
%
% stimo il trend 
A = vander(0:N-1); A=A(:,end-1:end);
b = log(y');
[Q, R] = qr(A);
xtrend = R \ (Q'*b);
figure(3); plot(0:N-1,exp(xtrend(2))*exp(xtrend(1)*(0:N-1)),'g-');  pause,
exp(xtrend(2))
xtrend(1)
pause,
%
% elimino il trend
yflutt =  y - exp(xtrend(2))*exp(xtrend(1)*(0:N-1));
figure(4); plot(0:N-1,yflutt,'r-'); pause,
%
figure(44); clf, hold on,
disp(['a = ' mat2str(a,4)]);
for na_est=3:na+5
  [U, S, V] = svd(hankel([zeros(1,na_est-1), yflutt(1:N-na_est+1)], yflutt(N-na_est+1:N-1)));
  figure(44); semilogy(log10(diag(S)),'bo'); diag(S);
  pause
  figure(44); semilogy(log10(diag(S)),'go'); 
  [Q, R] = qr(hankel([zeros(1,na_est-1), yflutt(1:N-na_est+1)], yflutt(N-na_est+1:N-1)));
  b = yflutt';
  xflutt = R \ (Q'*b);
  a_est = [1 -xflutt(end:-1:1)'];
  disp(['a_est = ' mat2str(a_est,4)]);
end;
%
% errore di predizione a un passo:
pe = zeros(1,length(y));
y_est = yflutt;
for i=na_est:length(y)
  y_est(i) = 1/a_est(1) * sum( -a_est(2:end).*yflutt(i-1:-1:i-na_est+1) );
  pe(i) = yflutt(i) - y_est(i);
end;
y_est = y_est + exp(xtrend(2))*exp(xtrend(1)*(0:N-1));
figure(3); plot(0:N-1,y_est,'b-'); pause,
figure(33); plot(0:N-1,pe,'g-'); pause,
disp(['varianza dell''errore di predizione a un passo: ' num2str(var(pe))]);
%
% vediamo se ho "sbiancato" la serie; applico il periodogramma cumulato di Bartlett ad "yflutt" ed a
% "pe" :
[Pff, f] = psd(yflutt);
%figure(4); clf, hold on, plot(f, 10*log10(Pff), 'r-');
est_mean_f = mean(yflutt);
est_var_f = sum((yflutt - est_mean_f).^2)/(N-1);
CPss_yf = zeros(length(f), 1);
for i=1:length(f)
  CPss_yf(i) = sum(Pff(1:i)) / (N * est_var_f);
end;
figure(5); clf, hold on, plot(f, CPss_yf, 'r-');
[Pee, f] = psd(pe);
%figure(4); plot(f, 10*log10(Pee), 'b-'); print -depsc LEZ_10_figure_4.epsc; 
est_mean_e = mean(pe);
est_var_e = sum((pe - est_mean_e).^2)/(N-1);
CPss_pe = zeros(length(f), 1);
for i=1:length(f)
  CPss_pe(i) = sum(Pee(1:i)) / (N * est_var_e);
end;
figure(5); plot(f, CPss_pe, 'b-');
% per confronto riporto il periodogramma cumulato del rumore bianco:
var_s = 6;
s = randn(1,N) * sqrt(var_s);
[Pss, f] = psd(s);
est_mean_s = mean(s);
est_var_s = sum((s - est_mean_s).^2)/(N-1);
CPss_wn = zeros(length(f), 1);
for i=1:length(f)
  CPss_wn(i) = sum(Pss(1:i)) / (N * est_var_s);
end;
figure(5); plot(f, CPss_wn, 'g-');


