%
clear;

%... read the data

sp500=importdata('homework03_SP500.txt');
infla=importdata('homework03_inflation.txt');
yield=importdata('homework03_yield.txt');
div=2.00;  %  https://www.multpl.com/s-p-500-dividend-yield
sp500 = sp500/sp500(1); spdiv(1)=1; dinf(1)=1; 
for i=2:length(sp500); 
    sprise=(sp500(i)-sp500(i-1))/sp500(i-1);
    spdiv(i)=spdiv(i-1)*(1+sprise)*(1+0.01*div/12);
    dinf(i)=dinf(i-1)*(1+infla(i)*0.01/12);
end

%  task 1 plot
lw=2;  figure(1);
p1=semilogy(sp500,'Linewidth',lw);

set(gca,'FontSize',25);  xlabel('months'); ylabel('index');
title('SP 500  etc');  set(gcf,'Position', [ 100 100 1000 1000 ]);
hold on;
p2=semilogy(spdiv,'Linewidth',lw); hold on;
p3=semilogy(dinf,'Linewidth',lw); hold on;


%...... plot linear yield and compound yield
b1(1)=1; b2(1)=1;  for i=2:length(sp500); 
    b1(i)=b1(i-1)+0.01*yield(i)/12;
    b2(i)=b2(i-1)*(1+0.01*yield(i)/12);
end
p4=semilogy(b1,'Linewidth',lw);  hold on; p5=semilogy(b2,'Linewidth',lw); hold off;
%   can only add legends at the end for multiple lines
legend([p1,p2,p3,p4,p5],["SP500","SP500/div","Inflation","Yield","Yield/compunded"],  ...
    'Location','Northwest');  hold on;
hold off; return;

%  task 2: (not) fit and add to plot
%  y=a*exp(b*x)  -->  fit log(y)=log(a) + b*x > a=10^a0; b=a1
%N=length(sp500); x=1:N; y=log10(sp500);  [c0,c1]=fit(x,y);  
%A=10.0^c0; B=c1;
%fprintf('SP500   fit:  A=%12.3f  B=%12.6f\n',A,B*12);
%y=A*10 .^ (B*x);  semilogy(y);


%... task 3:  only the histogram for different intervals 5y, 10y, 20y
%    
years=10;  clear h;  
k1=0; k2=0; rr1=0; rr2=0;  kk=0; 
for m1=1:860-12*years;  kk=kk+1; 
    [z1,z2]=invest(spdiv,m1,m1+12*years);  
    % yearly rate of return
    rate=exp( log(z1/(12*years))/years)-1;  rate=rate*100;  %  percent
    h(kk)=rate;
    zz=z1/z2; fprintf('  %12.2f  rate=%12.2f  zz=%12.5f\n',z1,rate,zz);
    rr1=rr1+z1;  rr2=rr2+z2;
    if zz>1; k1=k1+1; else; k2=k2+1; end
end
fprintf('  zz>0:  %d  zz<1:  %d  rr1=%f  rr2=%f\n',k1,k2,rr1/kk,rr2/kk);
figure(2);  histogram(h,101);  grid on;


%  task 2 function
function [r1,r2]=invest(s,m1,m2)
    r1=0; r2=0;
    %....  exp growth w/o dips
    %  ym=y0*10^am -->  a=log10(ym/y0)/m
    a=log10(s(m2)/s(m1))/(m2-m1+1);
    for m=m1:m2
        %.... ivest one dollar every month, accumate shares
        r1=r1+1.0/s(m);
        s2=s(m1)*10^(a*(m-m1));
        r2=r2+1.0/s2;
    end
    %... at the end we 'seel' the shares at the price of that month
    r1=r1*s(m2);  r2=r2*s(m2);
end



function [a0,a1]=fit(X,Y)
    N=length(X); xi=0; xi2=0; yi=0; xiyi=0;
    for i=1:N
        xi=xi+X(i);
        yi=yi+Y(i);
        xi2=xi2+X(i)*X(i);
        xiyi=xiyi+X(i)*Y(i);
    end
    a1=(N*xiyi-xi*yi)/(N*xi2-xi*xi);
    a0=(yi/N) - a1*xi/N;
end


