% Linear Prediction:given y(1),y(2),...,y(n-1), want to predict y(n)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%load sunspot
%load dow_jones
%Generating Articifical Data
n=(1:0.5:25)';
y=3*sin(2*pi*0.52*n)+2*sin(2*pi*0.40*n);
noise=input('Enter Residual Level 0 ~ 1');
y=y+noise*normrnd(0,1,length(y),1);%add residuals into the original data if necessary

[mm nn]=size(y);
if mm==1
    y=y';%make data a vertical vector
end
N=length(y);
fprintf('¸ê®Æªø«×¬° %7.0f \n',N);
m=input('¿é¤J¦ô­p¥Îªº¸ê®Æªø«×¡G');
p=input('Enter the order p for all approaches: ');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   
%   moving average-like approach
%   Model:  y(n)=(y(n-1)+y(n-2)+...+y(n-p))/p
%
A=toeplitz(y(p:m-1),y(p:-1:1));
if p>1
    y_hat=mean(A')';
else
    y_hat=A;
end
error1=norm(y(p+1:m)-y_hat)/norm(y(p+1:m))
%plot([y(p+1:m) y_hat])
%plot real data vs. fitted data
subplot(2,2,1),plot([y(p+1:m) y_hat])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   moving average-like approach with weightings
%   Model:  y(n)=(3/6)y(n-1)+(2/6)y(n-2)+(1/6)y(n-3)
%
A=toeplitz(y(3:m-1),y(3:-1:1));
y_hat2=A*[3/6 2/6 1/6]';
error2=norm(y(4:m)-y_hat2)/norm(y(4:m))
subplot(2,2,2),plot([y(4:m) y_hat2])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   moving average-like approach with degrading weightings
%   Model:  y(n)=alpha*y(n-1)+alpha(1-alpha)y(n-1)+alpha(1-alpha)^2y(n-1)+.....
%
alpha=0.99;
coef=alpha*ones(p,1).*(1-alpha).^(0:p-1)';
A=toeplitz(y(p:m-1),y(p:-1:1));
y_hat3=A*coef;
error3=norm(y(p+1:m)-y_hat3)/norm(y(p+1:m))
subplot(2,2,3),plot([y(p+1:m) y_hat3])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Model:   Least square error method
%
yy=y(p+1:m);
X=Toeplitz(y(p:m-1),y(p:-1:1));
a=pinv(X)*yy;
y_hat4=X*a;
error4=norm(y(p+1:m)-y_hat4)/norm(y(p+1:m))
subplot(2,2,4),plot([y(p+1:m) y_hat4])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Model:   near null space method
XX=[yy -X];
[m n]=size(XX);
%To generate an orthonormal set of basis for the Null(XX'*XX)
%(1) use eigenvectors
[V,S]=eig(XX'*XX);
diag(S)
s=input('Enter the number of zero eigenvalues :');%how many non-zero eigenvalues        
V_svd=V(:,1:s);
%%%%%%%%%%%%%%%%%%%%%
%USE SVD
%[U,S,V]=svd(XX'*XX);  
%diag(S)
%s=input('Enter the number of non-zero eigenvalues :');%how many non-zero eigenvalues        
%V_svd=V(:,s+1:n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[error_pre,error5,y_svd]=near_null(y,V_svd,m,p);
error5