clear
format compact
start=1;  % Starting value in first year
rate=.07; % Annual rate of return
Noise=.2;
NumTrials=30;
disp(['TrueRateOfReturn = ' num2str(rate)])
x=1:66; % time in years
y=start.*(1+rate).^x;
for trial=1:NumTrials,
    % Add proportional random noise to y data
    yy=y+Noise.*y.*randn(size(x));
    %  Method 1: fit noisy log y transfomed data to a straight line
    figure(1)
    [coeff,trR2]=plotit(x,log(abs(yy)),1);
    start1(trial)=exp(coeff(2));
    rate1(trial)=exp(coeff(1))-1;
    xlabel('Time, years')
    ylabel('Natural log of value')
    title('Straight-line fit to log(value) transformed data')
    % Method 2: use iterative line-linear fit to noisy x,y data
    figure(2)
    [rate2(trial),start2(trial),FittingError,fsR2]=fitshape1(x,yy,0);
    xlabel('Time, years')
    ylabel('Value')
    title('Iterative fit to y=start*(1+rate)^x to noisy x,y data')
end
meanstart1=mean(start1);stdstart1=std(start1);rsdstart1=stdstart1./meanstart1;
meanstart2=mean(start2);stdstart2=std(start2);rsdstart2=stdstart2./meanstart2;
meanrate1=mean(rate1);stdrate1=std(rate1);rsdrate1=100.*(stdrate1./meanrate1);
meanrate2=mean(rate2);stdrate2=std(rate2);rsdrate2=100.*(stdrate2./meanrate2);
range1=max(rate1)-min(rate1)
range2=max(rate2)-min(rate2)
text(min(x),max(y)-.05.*(max(y)-min(y)),['   Rate of return = ' num2str(meanrate2)] );
disp('                           Rate          % RSD    high-low range')
disp(['Coordinate transformation: ' num2str([meanrate1 rsdrate1 range1])])
disp(['Iterative curve fitting:   ' num2str([meanrate2 rsdrate2 range2])])
disp(' ')