'Figure 1'
[Date Temp] = textread('McGuireAFB.dat', '%s%f');
Time = 1955+1/2 + (1:length(Temp))'/365.25;
figure(1);
x = Time;
y = Temp;
plot(x,y,'b.');
lsline;
xlabel('Year');
ylabel('Temperature (F)');
title('Time Series');
xlim([1955 2010]);
n = length(x);
xbar = mean(x);
ybar = mean(y);
sx = std(x);
sy = std(y);
zx = (x-xbar)/sx;
zy = (y-ybar)/sy;
r = sum(zx.*zy)/(n-1);
b1 = r*sy/sx;
b0 = ybar - b1*xbar;
       sprintf('y = %5.3f + %5.3f x / 100', b0, 100*b1)
legend(sprintf('Temp = %5.3f + %5.3f t / 100', b0, 100*b1),'Location','NorthWest');
bigger;

'Figure 2'
window = ones(365,1)/365;
Tempw = conv(Temp,window,'valid');
Time = 1955+1/2 + (1:length(Tempw))'/365.25;
x = Time;
y = Tempw;
figure(2);
plot(x,y,'b.');
lsline;
xlabel(sprintf('Year (%4d + t)',round(maxx)));
ylabel('Temperature (F)');
title('Time Series');
xlim([1955 2010]);
n = length(x);
maxx = max(x);
x = x-max(x);
xbar = mean(x);
ybar = mean(y);
sx = std(x);
sy = std(y);
zx = (x-xbar)/sx;
zy = (y-ybar)/sy;
r = sum(zx.*zy)/(n-1);
b1 = r*sy/sx;
b0 = ybar - b1*xbar;
R2 = r^2;
mystr1 = sprintf('Temp = %5.3f + %5.3f t / 100', b0, 100*b1);
mystr2 = sprintf('R^2 = %f',R2)
legend(mystr1, mystr2, 'Location','NorthWest');
bigger;

'Figure 3'
figure(3);
x = Time(1:365:end);
y = Tempw(1:365:end);
plot(x,y,'r-');
hold on;
plot(x,y,'r+');
hold off;
xlabel(sprintf('Year (%4d + t)',round(maxx)));
ylabel('Temperature (F)');
title('Time Series');
lsline;
n = length(x);
maxx = max(x);
x = x-max(x);
xbar = mean(x);
ybar = mean(y);
sx = std(x);
sy = std(y);
zx = (x-xbar)/sx;
zy = (y-ybar)/sy;
r = sum(zx.*zy)/(n-1);
b1 = r*sy/sx;
b0 = ybar - b1*xbar;
R2 = r^2;
mystr1 = sprintf('Temp = %5.3f + %5.3f t / 100', b0, 100*b1);
mystr2 = sprintf('R^2 = %f',R2)
legend(mystr1, mystr2, 'Location','NorthWest');
bigger;

'Figure 4'
[yr mo day Temp] = textread('McGuireAFB.dat', '%4s%2s%2s%f');
figure(4);
boxplot(Temp,mo);
xlabel('Month');
ylabel('Temperature (F)');
title('Figure 3.7');
bigger;

'Figure 5'
[yr mo day Temp] = textread('McGuireAFB.dat', '%4d%2d%2d%f');
moday = (mo-1)*30 + day;
figure(5);
plot(moday,Temp,'k.');
xlabel('Day of the Year');
ylabel('Temperature (F)');
title('Figure 3.8');
bigger;

'Figure 6'
mo   = mo(1:365);
day  = day(1:365);
Temp = Temp(1:365);
moday = (mo-1)*30 + day;
figure(6);
plot(moday,Temp,'k.');
%
window = ones(20,1)/20;
Tempw = conv(Temp,window,'same');
hold on;
plot(moday,Tempw,'r-');
hold off;
%
window = ones(80,1)/80;
Tempw = conv(Temp,window,'same');
hold on;
plot(moday,Tempw,'b-');
hold off;
%
xlabel('Day of the Year');
ylabel('Temperature (F)');
title('Figure 3.8');
bigger;

'Figure 7'
figure(7);
[Date Temp] = textread('McGuireAFB.dat', '%s%f');
window = ones(365,1)/365;
Tempw = conv(Temp,window,'valid');
Time = 1955+1/2 + (1:length(Tempw))'/365.25;
x = Time(1:365:end);
y = Tempw(1:365:end);
n = length(x);
maxx = max(x);
x = x-max(x);
xbar = mean(x);
ybar = mean(y);
sx = std(x);
sy = std(y);
zx = (x-xbar)/sx;
zy = (y-ybar)/sy;
r = sum(zx.*zy)/(n-1);
b1 = r*sy/sx;
b0 = ybar - b1*xbar;
       sprintf('y = %5.3f + %5.3f x / 100', b0, 100*b1)
yhat = b0 + b1*x;
residuals = y-yhat;
plot(yhat,residuals,'r+');
xlabel('Predicted Temp (F)');
ylabel('Residual (F)');
title('Residuals for Figure 3');
%legend(sprintf('Temp = %5.3f + %5.3f (%4d+t) / 100', b0, 100*b1, round(maxx)),'Location','NorthWest');
bigger;

'Figure 8'
figure(8);
[Date Temp] = textread('McGuireAFB.dat', '%s%f');
window = ones(365,1)/365;
Tempw = conv(Temp,window,'valid');
Time = 1955+1/2 + (1:length(Tempw))'/365.25;
x = Time(1:365:end);
y = Tempw(1:365:end);
n = length(x);
maxx = max(x);
x = x-max(x);
xbar = mean(x);
ybar = mean(y);
sx = std(x);
sy = std(y);
zx = (x-xbar)/sx;
zy = (y-ybar)/sy;
r = sum(zx.*zy)/(n-1);
b1 = r*sy/sx;
b0 = ybar - b1*xbar;
       sprintf('y = %5.3f + %5.3f x / 100', b0, 100*b1)
yhat = b0 + b1*x;
residuals = y-yhat;
plot(x,residuals,'r+');
xlabel('Year');
ylabel('Residual (F)');
title('Residuals for Figure 3');
%legend(sprintf('Temp = %5.3f + %5.3f (%4d+t) / 100', b0, 100*b1, round(maxx)),'Location','NorthWest');
bigger;

'Figures 9 and 10'
[Date Temp] = textread('McGuireAFB.dat', '%s%f');
window = ones(365,1)/365;
Tempw = conv(Temp,window,'valid');
Time = 1955+1/2 + (1:length(Tempw))'/365.25;
x = Time(1:365:end);
y = Tempw(1:365:end);
n = length(x);
maxx = max(x);
x = x-max(x);
xbar = mean(x);
ybar = mean(y);
sx = std(x);
sy = std(y);
zx = (x-xbar)/sx;
zy = (y-ybar)/sy;
r = sum(zx.*zy)/(n-1);
b1 = r*sy/sx;
b0 = ybar - b1*xbar;
b0b1orig = sprintf('y = %5.3f + %5.3f x / 100', b0, 100*b1)
yhat = b0 + b1*x;
residuals = y-yhat;
for iter=1:1000
    indices = random('Discrete Uniform',n,n,1);
    y = yhat + residuals(indices);
    ybar = mean(y);
    sy = std(y);
    zy = (y-ybar)/sy;
    r = sum(zx.*zy)/(n-1);
    betahat1(iter) = r*sy/sx;
    betahat0(iter) = ybar - betahat1(iter)*xbar;
end
figure(9);
histogram(betahat0);
xlabel('b_0 in degrees (F)');
ylabel('Bootstrap frequency');
title('Histogram of b_0 values');
bigger;
figure(10);
histogram(100*betahat1);
xlabel('b_1 in degrees (F) per century');
ylabel('Bootstrap frequency');
title('Histogram of b_1 values');
bigger;

'Figure 11'
[Date Temp] = textread('McGuireAFB.dat', '%s%f');
Time = 1955 + (1:length(Temp))'/365.25;
minTime = min(Time);
Time = Time-minTime;
figure(11);
n = length(Time);
X = [ones(n,1) Time cos(2*pi*Time) sin(2*pi*Time)];
y = Temp;
plot(minTime+Time,Temp,'b.');
b = X\y;
ybar = X*b;
hold on;
plot(minTime+Time,ybar,'r-');
hold off;
xlabel('Year (1955 + t)');
ylabel('Temperature (F)');
title('Time Series');
R2 = 1 - sum((y-X*b).^2)/sum((y-mean(y)).^2);
mystr1 = sprintf('Temp = %5.3f + %5.3f t / 100 + %5.3f cos(2 \\pi t) + %5.3f sin(2 \\pi t)', ...
	b(1), 100*b(2), b(3), b(4));
mystr2 = sprintf('R^2 = %f',R2)
legend(mystr1, mystr2, 'Location','NorthWest');
bigger;

'Figure 12'
figure(12);
residuals = y-X*b;
plot(minTime+Time,residuals,'b.');
xlabel('Year');
ylabel('Residual (degrees F)');
title('Residuals');
bigger;

'Figures 13 and 14'
n = length(y);
yhat = X*b;
residuals = y-yhat;
for iter=1:1000
    indices = random('Discrete Uniform',n,n,1);
    ynew = yhat + residuals(indices);
    b = X\ynew;
    coldest(iter) = 1+365.25*atan(b(4)/b(3))/(2*pi);
    b2(iter) = 100*b(2);
end
figure(13);
histogram(b2);
xlabel('Slope in Degrees (F) per Century');
ylabel('Bootstrap frequency');
title('Histogram for Warming/Cooling Trend');
bigger;
figure(14);
histogram(coldest);
xlabel('Coldest day of the year');
ylabel('Bootstrap frequency');
title('Histogram for Coldest Day');
bigger;