
Y=[9,15,30,40,56,66,84,102,131,159,173,186,190,221,248,278,330,354,382,461,481,526,587,608,697,781,896,999,1124,1212,1385,1715,2055,2429,2764,3332,4288,5373,6789,8564,10298]';
Y=[Y(1);diff(Y)];
T=[1:length(Y)]';
logY=log(Y);
A = [ones(length(Y),1) T];
a = A\logY

w = logY - (A*a);
stddev = std(w);
upper = exp(a(1)+a(2)*T+2*stddev);
middle = exp(a(1)+a(2)*T);
lower = exp(a(1)+a(2)*T-2*stddev);

plot(T,upper,':');
hold on
plot(T,lower,':');
hold on
plot(T,middle);
hold on
plot(T,Y)
%Problem 3
T1=[1:50]';
w0 = logY - (A*a);
stddev0 = std(w0);
upper0 = exp(a(1)+a(2)*T1+2*stddev0);
middle0 = exp(a(1)+a(2)*T1);
lower0 = exp(a(1)+a(2)*T1-2*stddev0);
figure;
plot(T1,cumsum(upper0),':');
hold on
plot(T1,cumsum(lower0),':');
hold on
plot(T1,cumsum(middle0));
hold on
plot(T,cumsum(Y));

DL=[0;0];
for i=1:length(Y)
    temp=[-exp(a(1)+a(2)*i)+Y(i);-i*exp(a(1)+a(2)*i)+i*Y(i)];
    DL=DL+temp;
end
DL=DL

HL=[0 0;0 0];
for i=1:length(Y)
    temp2=[-exp(a(1)+a(2)*i) -i*exp(a(1)+a(2)*i);-i*exp(a(1)+a(2)*i) -(i^2)*exp(a(1)+a(2)*i)];
    HL=HL+temp2;
end
HL=HL

a0 = [2.12;0.07];

for i=1:50
    DL=[0;0];
    HL=[0 0;0 0];
    for i=1:length(Y)
        temp=[-exp(a0(1)+a0(2)*i)+Y(i);-i*exp(a0(1)+a0(2)*i)+i*Y(i)];
        DL=DL+temp;
    end
    for i=1:length(Y)
        temp2=[-exp(a0(1)+a0(2)*i) -i*exp(a0(1)+a0(2)*i);-i*exp(a0(1)+a0(2)*i) -(i^2)*exp(a0(1)+a0(2)*i)];
        HL=HL+temp2;
    end
    a0=a0-HL\DL;
end
a0=a0

upper1 = exp(a0(1)+a0(2)*T)+2*sqrt(exp(a0(1)+a0(2)*T)); 
lower1 = exp(a0(1)+a0(2)*T)-2*sqrt(exp(a0(1)+a0(2)*T));
middle1 = exp(a0(1)+a0(2)*T);

figure;
plot(T,upper1,':');
hold on
plot(T,lower1,':');
hold on
plot(T,middle1);
hold on
plot(T,Y);

figure;
plot(T,cumsum(upper1),':');
hold on
plot(T,cumsum(lower1),':');
hold on
plot(T,cumsum(middle1));
hold on
plot(T,cumsum(Y));

T=[1:200];
Y2=[9,15,30,40,56,66,84,102,131,159,173,186,190,221,248,278,330,354,382,461,481,526,587,608,697,781,896,999,1124,1212,1385,1715,2055,2429,2764,3332,4288,5373,6789,8564,10298]';
T2=[1:length(Y2)];

upper = exp(a(1)+a(2)*T+2*stddev);
middle = exp(a(1)+a(2)*T);
lower = exp(a(1)+a(2)*T-2*stddev);

upper1 = exp(a0(1)+a0(2)*T)+2*sqrt(exp(a0(1)+a0(2)*T)); 
lower1 = exp(a0(1)+a0(2)*T)-2*sqrt(exp(a0(1)+a0(2)*T));
middle1 = exp(a0(1)+a0(2)*T);

figure;
semilogy(T,cumsum(upper1),':');
hold on
semilogy(T,cumsum(lower1),':');
hold on
semilogy(T,cumsum(middle1));
hold on
semilogy(T2,Y2);

hold on
semilogy(T,cumsum(upper),':');
hold on
semilogy(T,cumsum(lower),':');
hold on
semilogy(T,cumsum(middle));
hold on
semilogy(T,cumsum(Y));