%PROGRAM NAME: pete_stocks.m
%This code uses data from test_of_daily_suite.m
nw=10; %Moving Average Window Length
nx=length(xpr);
xMA=zeros(nx,1);
for k=nw:nx
    xMA(k)=mean(xpr(k-nw+1:k));
end
%MA FILTER:
nz=[1,zeros(1,nw-1),-1];
dz=[1,-1,zeros(1,nw-1)];
MAF=(1/nw)*tf(nz,dz,1);
t=1:nx; t=t';
%Use 1st order LPF to smooth:
wBW=.05*pi; %Use PSD(x) to choose
s=tf('s');
k0=2; %Set k0=2 for 1st order filter & k0=3 for 2nd order filter
if k0==3;
    wBW=3*wBW;
end
kk=0;
if k0==2
    kk=1;
end
LPF=(1/(1+s/wBW)^(k0-1));
LPFz=c2d(LPF,1,'imp');
[num,den]=tfdata(LPFz,'v');
nsum=sum(num); dsum=sum(den); %TF num and den coefficients
num=num*(dsum/nsum); %Adjust to unity static gain
xLPF=zeros(nx,1);
for k=k0:nx
    xLPF(k)=-den(2:k0)*xLPF(k-1:-1:k-k0+1)+num(2-kk:k0-kk)*xpr(k-1+kk:-1:k-k0+1+kk);
end
figure(1)
plot(t,[xpr,xMA,xLPF])
xlabel('Time (day #)')
title('GOOGL Close Price,and MA/LPF Estimates')
grid
legend('data','MA Estimate','LPF Estimate')
%--------------------------------
%MA and LPR Bode Plots
figure(2)
bode(MAF,{10^-1,10})
hold on
bode(LPFz,{10^-2,10},'r')
title('MA and 1ST Order Filter Bode Plots')
legend('MAF','LPF')
grid
%--------------------------------
%Residuals:
yMA=xpr-xMA; yLPF=xpr-xLPF;
figure(3)
plot(t,[yMA,yLPF])
xlabel('Time (day #)')
title('MA/LPF Estimate Residuals')
grid
legend('MA Estimate Residual','LPF Estimate Residual')
%====================================================
%====================================================
%1-Day Ahead Prediction using MA and LPF
xhatLPF=zeros(nx,1);
xhatMA=zeros(nx,1);
nw=10;
for k=nw+2:nx-1    
   x1=xpr(k-nw-1:k-1);
   x2=xpr(k-nw:k);
   C=cov(x1,x2);
   mLPF=C(1,2)/C(1,1);
   bLPF=mean(x2)-mLPF*mean(x1);
   xhatLPF(k+1)=mLPF*xpr(k)+bLPF;
   xhatMA(k+1)=xMA(k);
end
errLPF=xpr-xhatLPF;
errMA=xpr-xhatMA;
zvec=1390:1400; nzvec=length(zvec);
errLPF(zvec)=zeros(nzvec,1); errMA(zvec)=zeros(nzvec,1);
figure(4)
plot(t(2:nx),[errMA(2:nx),errLPF(2:nx)])
title('1-Day Ahead MA and LPF Prediction Error')
legend('MA error','LPF error')
grid
rmseMA=sqrt(sum(errMA(20:nx).^2))
rmseLPF=sqrt(sum(errLPF(20:nx).^2))
%========================================
%========================================
% Residual LINEAR PREDICTION
y=errLPF;
yhat=zeros(nx,1);
nw=7; % ***** Correlation Window Width ******
m=zeros(nx,1);
for k=nw+2:nx-1    
   y1=y(k-nw-1:k-1);
   y2=y(k-nw:k);
   C=cov(y1,y2);
   m(k)=C(1,2)/C(1,1);
   yhat(k+1)=m(k)*y(k);
end
erry=y-yhat;
figure(5)
plot(t,[y,yhat,erry])
xlabel('Time (day #)')
title('Residual 1-Day Ahead Prediction')
legend('LPF Residual','Residual Prediction','Error')
grid
%FINAL PREDICTION:
xhat=xhatLPF+yhat;
figure(6)
plot(t,[xpr,xhat])
title('Stock Price and 1-Day-Ahead Predicted Price')
legend('Close Price','LPF Predicted Close Price','Location','SouthEast')
grid
figure(5)
figure(4)
figure(3)
figure(2)
figure(1)
%Suppose that at day 1450 I have 100 shares valued at $670/share
n0=1450;
profit=zeros(nx-n0+1,1);
N=nx-n0;
j=1;
for k=n0:nx-1
    sell=0;
    if xhat(k+1)-xpr(k) >= 5
        sell=(xpr(k+1)-xpr(k))*10; %Sell 10 shares @ end of day k+1
    end
    profit(j+1)=profit(j)+sell;
    j=j+1;
end
figure(11)
plot(profit)
title('Profit From Selling 10 Shares if xhat(k+1)-x(k) > $10')
grid

        
