load SandPfrom1950

% T=TTT(1:75)'; % Time in years
% V=YYY(1:75)'; % S&P value

 T=x';
V=y';

% Method 1: iterative exponential fit to time,value data, using the
% fitshape1 function with the expression for compound interest in line 106.
figure(1)
[rate1,start1,FittingError,fsR2]=fitshape1(T,V,0);
xlabel('Time, years')
ylabel('S&P Index')
title('Iterative fit to y=start*(1+rate)^x to time,value data')
text(min(T),max(V)-.05.*(max(V)-min(V)),['   Rate of return = ' num2str(rate1)] );
grid
hold on;
plot(TTT,start1.*(1+rate1).^TTT,'b-')
hold off

% Method 2: fit log(value) transfomed data to a straight line
figure(2)
[coeff,trR2]=plotit(T,log(V),1);
start2=exp(coeff(2));
rate2=exp(coeff(1))-1;
xlabel('Time, years')
ylabel('Natural log of S&P Index')
title('Straight-line fit to time,log(value) transformed data')
grid 
disp(['Iterative curve fitting:   ' num2str(rate1)])
disp(['Coordinate transformation: ' num2str(rate2)])

function [a,b,FittingError,R2]=fitshape1(x,y,start)
% Simplified peak fitting example using fminsearch function. Fits a model
% shape defined the the shape function (line 60) to the data in x,y, with
% the first-guess values defined by vector 'start', in the form [position1
% position2, etc]. The shape must mave one iterated variable (e.g.
% position). Number of peaks in the model equals the length of 'start'.
% Change the last line to the expression for your shape. 
% Version 1, T. C. O'Haver, toh@umd.edu, 2016
%
% Example 1: Compound Interest
% NOTE: Change the last line to the expression for compound interest. 
% x=1:66; % time in years
% rate=.07; % Annual rate of return
% Noise=.5;
% y=start.*(1+rate).^x; %  Compound Interest formula
% yy=y+Noise.*y.*randn(size(x)); % Add proportional noise
% rate=fitshape1(x,yy,0)
%
format compact
global B
% Call the fminsearch function with specified fitting function
options = optimset('TolX',.000001,'TolFun',.00001,'Display','off','MaxFunEvals',1000);
FitResults=fminsearch(@(lambda)(fitfunction(lambda,x,y)),start,options);
NumPeaks=round(length(start));

% Construct model from FitResults
for m=1:NumPeaks
  A(m,:)=shapefunction(x,FitResults(m));
end
model=B'*A;

% Compute percent fitting error 
FittingError=100*norm(y-model)./(sqrt(length(x))*max(y));
% Compute R2
SStot=sum((y-mean(y)).^2);
SSres=sum((y-model).^2);
R2=1-(SSres./SStot);

% Plot data as dots, fitted model as line
plot(x,y,'r.',x,model,'b-')
hold on;
for m=1:NumPeaks
    plot(x,A(m,:)*B(m),'g');
end
hold off
title('fitshape1')
for m=1:NumPeaks
  a(m)=FitResults(m);
end
b=B';
end
% -----------------------------------------------------

function err = fitfunction(lambda,t,y)
%  Fitting function for a signal consisting of overlapping peaks.
global B
A = zeros(length(t),round(length(lambda)));
for j = 1:length(lambda)
    A(:,j) = shapefunction(t,lambda(j))';
end
B = A\y';
z = A*B;
err = norm(z-y');
end
% ----------------------------------------------------

function g = shapefunction(x,a)
%  shapefunction(x,a)
%  x may be scalar, vector, or matrix, a is scalar
%  T. C. O'Haver, 2016
% COMPOUND INTEREST (a=interest rate)
% g=(1+a).^x;
%
g=(1+a).^x;
end