%Thompson6B (uses mvtcdf)
format compact

%Set temperature increment dx:

dx = .1 % for testing program  Takes 19 seconds on desktop
dx = .01 % for production  Takes 190 seconds on desktop, 609 sec. on laptop.
dx 

% set quantiles to evaluate: 
pvalue = [.025 .25 .5 .75 .975];

% Read 5 cores from Climatic Change 2003 Fig 5 table file:
% source: bprc.osu.edu/Icecore/Climatic-change-2003-Fig5-table.XLS, 
% converted to .TXT with headers deleted and NaN for missing.
% Guliya col. 2
% Dunde col. 4
% Dasuopu col. 6
% (Quelccaya Summit missing)
% Sajama col. 9
% Huascaran Core 2 col. 13
% reverse order (1990s first)
% Guliya, Dunde, Dasuopu use 1-decades (1991-2000 etc)
% Sajama, Huascaran2 use 0-decades (1990-1999 etc) (Quelccaya as well)
% Assume 1-decades throughout. 
x = dlmread('Climatic-change-2003-Fig5-table.txt');
size(x) % ~ 100 x 16
d18O = NaN(100,6);
d18O(:,1) = x(:,2);
d18O(:,2) = x(:,4);
d18O(:,3) = x(:,6); 
d18O(:,5) = x(:,9);  % Quelccaya Summit goes in column 4
d18O(:,6) = x(:,13);
dec = (1991:-10:1001)';  % readings are in reverse calendar order.  
% Read Quelccaya Summit Core data from 
% ftp://ftp.ncdc.noaa.gov/pub/data/paleo/icecore/trop/quelccaya/q83summ.txt
% Headers and footers deleted in file q83summA.txt  
% Col. 1 - year from 1984 down to 744
% col. 6 - d18O (o/oo)
% 8 columns in all - Matlab reads 9, but last is zeros.
x = dlmread('q83summA.txt');
size(x) % ~ 1241 x 9
% Previous program found that 0-decades of Quelccaya Summit 
% core fit Fig. 6 exactly, while 1-decades and Quelccaya core 1 
% did not.  
% 1979 is row 6, 1000 is row 985:
Q = x(6:985,6);
% take decadal averages (0-decades) and insert in col. 4 of d18O:
Q = reshape(Q,10,98);
Q = mean(Q)';
d18O(3:100,4) = Q;
% insert 1990s and 1980s read visually from Fig. 6:
d18O(1,4) = -16.64;
d18O(2,4) = -17.67; 
dec = flipud(dec);  % reverse dec and d18O to forward calendar order.
d18O = flipud(d18O);

figure
plot(dec,d18O(:,1:3))
title ('Himalayan cores from Thompson et al CC 03')
xlabel('Bottom year of decade')
ylabel('d18O (o/oo)')
legend('Guliya', 'Dunde', 'Dasuopu')
grid on

% B&W version, in Dunde (#2), Guliya (#1), Dasuopu (#3) order 
%  from top to bottom.
figure
hold on
plot(dec,d18O(:,2),'k-', 'LineWidth', 1)
plot(dec,d18O(:,1),'k-', 'LineWidth', 2)
plot(dec,d18O(:,3),'k-', 'LineWidth', 2.5)
title ('Data on d18O for three Himalayan cores')
xlabel('Bottom year of decade')
ylabel('d18O (o/oo)')
legend('Dunde', 'Guliya', 'Dasuopu')
grid on

figure
plot(dec,d18O(:,4:6))
title ('Andean cores from Thompson et al CC 03')
xlabel('Bottom year of decade')
ylabel('d18O (o/oo)')
legend('Quelccaya Summit', 'Sajama', 'Huascaran 2')
grid on

% B&W version, in Sajama (#5), Quelccaya (#4), Husacaran (#6) order, 
%  roughly top to bottom
figure
hold on
plot(dec,d18O(:,5),'k-')
plot(dec,d18O(:,4),'k-', 'LineWidth',2)
plot(dec,d18O(:,6),'k-', 'LineWidth',2.5)
title ('Data on d18O for three Andean cores')
xlabel('Bottom year of decade')
ylabel('d18O (o/oo)')
legend('Sajama', 'Quelccaya Summit', 'Huascaran 2')
grid on

means = nanmean(d18O);
a18O = d18O - repmat(means,100,1);
figure
plot(dec,a18O(:, 1:3))
title ('Demeaned d18O for Himalayan cores')
xlabel('Bottom year of decade')
ylabel('d18O (o/oo), demeaned')
legend('Guliya', 'Dunde', 'Dasuopu')
grid on

figure
plot(dec,a18O(:,4:6))
title ('Demeaned d18O for Andean cores')
xlabel('Bottom year of decade')
ylabel('d18O (o/oo), demeaned')
legend('Quelccaya Summit', 'Sajama', 'Huascaran 2')
grid on

std = nanstd(a18O);
z = a18O./repmat(std,100,1);
figure
plot(dec,z(:,1:3))
title ('Z scores for Himalayan cores')
xlabel('Bottom year of decade')
ylabel('z-score of d18O')
legend('Guliya', 'Dunde', 'Dasuopu')
grid on

figure
plot(dec,z(:,4:6))
title ('Z scores for Andean cores')
xlabel('Bottom year of decade')
ylabel('z-score of d18O')
legend('Quelccaya Summit', 'Sajama', 'Huascaran 2')
grid on

hcz = mean(z(:,1:3)')';  % mean skips NaNs 
acz = mean(z(:,4:5)')';
lcz = (acz + hcz)/2;     % LCZ as constructed by Thompson et al fig 7c.
lcz4 = (acz + z(:,3))/2;
lcz4N = [NaN(98,1); lcz(99); lcz4(100)];

figure
plot(dec, hcz, 'b', dec, z(:,3), 'r')
title ('Himalayan Composite Z-Score')
xlabel('Bottom year of decade')
ylabel('Composite z-score')
legend ('HCZ (3 cores)', 'Dasuopu only')
grid on

figure
plot(dec, hcz, 'b', dec(99:100), [hcz(99); z(100,3)], 'r')
title ('Himalayan Composite Z-Score')
xlabel('Bottom year of decade')
ylabel('Composite z-score')
legend ('HCZ (3 cores)', 'Dasuopu only')
grid on

figure
plot(dec,acz)
title ('Andean Composite Z-Score')
xlabel('Bottom year of decade')
ylabel('Composite z-score')
legend ('ACZ (3 cores)')
grid on

figure
plot(dec,lcz, 'b', dec, lcz4, 'r')
title ('Composite Z-Scores')
xlabel('Bottom year of decade')
ylabel('Composite z-score')
legend ('6-core LCZ', '4-core LCZ')
grid on 

figure
plot(dec,lcz, 'b', dec, lcz4N, 'r')
title ('Low-Latitude Composite Z-Scores')
xlabel('Bottom year of decade')
ylabel('Composite z-score')
legend ('6-core LCZ', '4-core LCZ')
grid on

% B&W version, with x for LCZ4
figure
plot(dec,[lcz(1:99); lcz4(100)], 'k-',dec(100),lcz4(100),'k-o')
title ('Low-Latitude Composite Z-Scores')
xlabel('Bottom year of decade')
ylabel('Composite Z-score')
legend ('6-core LCZ', '4-core LCZ')
grid on

x = dlmread('crutem3vgl.txt');
% source www.cru.uea.uk/cru/data/temperature/crutem3vgl.txt,
% Variance adjusted land air temperature anomalies, 1961-90 = 0.  
% Global average, Downloaded 11/17/09
% Cite Brohan, P., J.J. Kennedy, I. Harris, S.F.B. Tett and P.D. Jones, 
%  2006: Uncertainty estimates in regional and global observed 
%  temperature changes: a new dataset from 1850. J. Geophysical 
%  Research 111, D12106, doi:10.1029/2005JD006548 -- Available as PDF
% 1/1850 - 9/2009, last 3 mo of 2009 filled with 0's, but has ann. avg.  
% Each year has two rows:
%   year, 12 months, annual avg.  (14 numbers)
%   year, % coverage, blank.  (13 numbers, but blank read as 0 for col. 14)
% 160 years in all (320 rows).
disp('Size of crutem3vgl.txt       ')
size(x)  % ~ 320 x 14
% isolate annual averages: 
x = reshape(x', 28, 160)';
size(x)
cruyr = x(:,1);
cruann = x(:, 14);
figure
plot(cruyr, cruann)
title ('CRUTEM3vGL, annual averages')
xlabel('Year')
ylabel('Temperature anomaly, deg.C, 1961-90 = 0')
grid on
disp('mean for 2001-2009')
cru2001dec = mean(cruann(152:160))
% isolate 1851-2000:
cruann = cruann(2:151);
% calculate 1-decadal averages:
crutem = reshape(cruann,10,15);
crutem = mean(crutem)';
crudec = 1851:10:2000;
figure
plot(crudec, crutem)
title ('CRUTEM3vGL, decadal averages through 1991-2000')
grid on
ylabel('Temperature anomalies, deg. C, 1961-90 = 0')
xlabel('Bottom year of decade')

disp('LCZ6 on CRUTEM3vGL, 1851s - 1981s (14 obs)')
[b6, Cov6, se, tstat, pb, R2, s26, e, r1, dw, pdw] = ...
    OLS(lcz(86:99), [ones(14,1) crutem(1:14)]);
[b6 se tstat pb]
R2
s26
Cov6(1,2)
r1
dw
pdw

disp('LCZ4 on CRUTEM3vGL, 1851s - 1991s (14 obs)')
[b4, Cov4, se, tstat, pb, R2, s24, e, r1, dw, pdw] = ...
    OLS(lcz4(86:100), [ones(15,1) crutem(1:15)]);
[b4 se tstat pb]
R2
s24 
Cov4(1,2)
r1
dw
pdw


disp('Start Calibration tic toc')
tic
xord = -5:dx:5;
[recon, reconcdf] = calib(lcz, b6, Cov6, s26, (14-2), xord);
toc
figure
plot(dec, recon);
title ('6-core LCZ on CRUTEM3vGL')
grid on

figure
plot(xord, reconcdf')
title ('Recon CDF for 6-core LCZ on CRUTEM3vGL.')
ylabel ('Cumulative Probability')
xlabel ('Temp anomaly, deg. C (1961-90 = 0)')
grid on

disp('quantile tic')
tic
reconq = quantile(pvalue, reconcdf, xord);
toc
% .02 seconds with n = 100, length(pvalue) = 5
figure
plot(dec,[recon reconq])
title ('6-core LCZ on CRUTEM3vGL.  Recon, 50% CI, 95% CI')
ylabel ('Temp anomaly, deg. C (1961-90 = 0)')
xlabel ('Bottom year of decade')
legend ('recon', '.025', '.25', '.5', '.75', '.975')
grid on


xord = -5:dx:5;
h = size(xord,2);
recon4 = NaN(100,1);
reconcdf4 = NaN(100,h);
% restrict LCZ4 recon to 1851-2000 by setting r1 = 86.  r1 = 1 to 
% compute full recon.  
r1 = 1; % doubles time unnecessarily
r1 = 86
[recon4(r1:100), reconcdf4(r1:100,:)] = calib(lcz4(r1:100), b4, Cov4, s24, (15-2), xord);
figure
plot(dec, recon4);
title ('4-core LCZ to CRUTEM3vGL')
grid on

figure
plot(xord, reconcdf4')
title ('Recon CDF for 4-core LCZ on CRUTEM3vGL.')
ylabel ('Cumulative Probability')
xlabel ('Temp anomaly, deg. C (1961-90 = 0)')
grid on

reconq4 = quantile(pvalue, reconcdf4, xord);
disp('size(reconq4))')
size(reconq4)

figure
plot(dec,[recon4 reconq4])
title ('4-core LCZ to CRUTEM3vGL.  Recon, 50% CI, 95% CI')
ylabel ('Temp anomaly, deg. C (1961-90 = 0)')
xlabel ('Bottom year of decade')
legend ('recon', '.025', '.25', '.5', '.75', '.975')
xlim([1850 2000])
grid on

figure
plot(dec(86:100), [crutem recon(86:100) reconq(86:100,3) recon4(86:100) reconq4(86:100,3)])
% recon4 ~ 100x 100!
title ('CRUTEM3vGL, 6-core recon, 4-core recon');
legend('CRUTEM3vGL', '6-core recon', '6-core median', '4-core recon', ...
    '4-core median')
ylabel ('Temp. anomaly, deg. C (1961-90 = 0)')
xlabel ('Bottom year of decade')
grid on

%B&W version
figure
hold on
plot(dec(86:100), crutem,'k-o', ... 
     dec(86:100), recon(86:100), 'k-s', ... 
     dec(86:100), reconq(86:100,3), 'k--', ... 
     dec(86:100), recon4(86:100),'k-^', ... 
     dec(86:100), reconq4(86:100,3), 'k-.')
title ('Instrumental vs. Reconstructed Temperature');
legend('Instrumental', '6-core classical', '6-core median', '4-core classical', ...
    '4-core median')
ylabel ('Temp. anomaly, deg. C (1961-90 = 0)')
xlabel ('Bottom year of decade')
grid on
box on



% eliminate all but last observations on recon4 etc, replace 
%  next to last with recon etc:
recon4N = recon4;
recon4N(1:98) = NaN(98,1);
recon4N(99) = recon(99);
reconq4N = reconq4;
reconq4N(1:98,:) = NaN(98,5);
reconq4N(99,:) = reconq(99,:);
figure
plot(dec,recon,'b',dec,recon4N,'r')
title ('Temperature reconstruction')
legend ('6-core recon', '4-core recon')
ylabel ('Temp. anomaly, deg. C (1961-90 = 0)')
xlabel ('Bottom year of decade')
grid on

figure
plot(dec, recon, 'b', dec, recon4N, 'r', dec, ...
    reconq(:,[1 2 4 5]), 'b', dec, reconq4N(:,[1 2 4 5]), 'r')
title ('Temperature reconstruction with 50% CI, 95% CI')
legend('6-core recon, 50% CI, 95% CI', '4-core recon, 50% CI, 95% CI')
ylabel ('Temp. anomaly, deg. C (1961-90 = 0)')
xlabel ('Bottom year of decade')
grid on

% B&W version
figure
plot(dec, [recon(1:99); recon4(100)], 'k',... 
     dec(100), recon4(100), 'k-o',...
     dec, [reconq(1:99,[1 2 4 5]); reconq4(100, [1 2 4 5])], 'k', ...
     dec(100), reconq4(100,[1 2 4 5]), 'k-o')
title ('Temperature reconstruction with 50% and 95% CI')
legend('6-core reconstruction', '4-core reconstruction')
ylabel ('Temp. anomaly, deg. C (1961-90 = 0)')
xlabel ('Bottom year of decade')
grid on

[xprime, x1, x2, x3, x4] = fieller(lcz(1:99), b6, Cov6, s26, (14-2), .50);
[xprime(100), x1(100), x2(100), x3(100), x4(100)] = ...
    fieller(lcz4(100), b4, Cov4, s24, (15-2), .50);
figure
title ('Posterior vs Fieller 50% Confidence Intervals');
x1 = max(x1,-5);
x2 = max(x1,-5);
x3 = max(x3,-5);
x4 = max(x4,-5);
x1 = min(x1,5);
x2 = min(x2,5);
x3 = min(x3,5);
x4 = min(x4,5);
hold on
plot(dec, [reconq(1:99,2); reconq4(100, 2)], 'k') %Must plot, then overwrite this 
  % for legend to work right.
fill([dec; flipud(dec)], [x1; flipud(x2)], [.8 .8 .8], 'EdgeColor', [.8 .8 .8])
fill([dec; flipud(dec)], [x3; flipud(x4)], [.8 .8 .8], 'EdgeColor', [.8 .8 .8])
plot(dec, [reconq(1:99,[ 2 4]); reconq4(100, [ 2 4 ])], 'k', ...
     dec(100), reconq4(100,[ 2 4]), 'k-o')
legend ('Posterior', 'Fieller');
hold off
grid on
box on

% Color version
figure
title ('Posterior vs Fieller 50% Confidence Intervals');
hold on
plot(dec, [reconq(1:99,2); reconq4(100, 2)], 'b') %Must plot, then overwrite this 
  % for legend to work right.
fill([dec; flipud(dec)], [x1; flipud(x2)], [1 .7 .7], 'EdgeColor', [1 .7 .7])
fill([dec; flipud(dec)], [x3; flipud(x4)], [1 .7 .7], 'EdgeColor', [1 .7 .7])
plot(dec, [reconq(1:99,[ 2 4]); reconq4(100, [ 2 4 ])], 'b', ...
     dec(100), reconq4(100,[ 2 4]), 'b-o')
legend ('Posterior', 'Fieller');
hold off
grid on
box on

[zprime, z1, z2, z3, z4] = fieller(lcz(1:99), b6, Cov6, s26, (14-2), .05);
[zprime(100), z1(100), z2(100), z3(100), z4(100)] = ...
    fieller(lcz4(100), b4, Cov4, s24, (15-2), .05);
figure
title ('Posterior vs Fieller 95% Confidence Regions')
ymax = 4.97;
z1 = max(z1,-ymax);
z2 = max(z2,-ymax);
z3 = max(z3,-ymax);
z4 = max(z4,-ymax);
z1 = min(z1,ymax);
z2 = min(z2,ymax);
z3 = min(z3,ymax);
z4 = min(z4,ymax);
hold on
plot(dec, [reconq(1:99,1); reconq4(100,1)], 'k')
fill([dec; flipud(dec)], [z1; flipud(z2)], [.8 .8 .8], 'EdgeColor', [.8 .8 .8])
fill([dec; flipud(dec)], [z3; flipud(z4)], [.8 .8 .8], 'EdgeColor', [.8 .8 .8])
plot(dec, [reconq(1:99,[1 5]); reconq4(100, [1 5])], 'k', ...
     dec(100), reconq4(100,[1 5]), 'k-o') 
legend ('Posterior', 'Fieller');
grid on
box on
hold off
xlim([998 2000])

% color version
figure
title ('Posterior vs Fieller 95% Confidence Regions')
hold on
plot(dec, [reconq(1:99,1); reconq4(100,1)], 'b')
fill([dec; flipud(dec)], [z1; flipud(z2)], [1 .7 .7], 'EdgeColor', [1 .7 .7])
fill([dec; flipud(dec)], [z3; flipud(z4)], [1 .7 .7], 'EdgeColor', [1 .7 .7])
plot(dec, [reconq(1:99,[1 5]); reconq4(100, [1 5])], 'b', ...
     dec(100), reconq4(100,[1 5]), 'b-o') 
legend ('Posterior', 'Fieller');
grid on
box on
hold off
xlim([998 2000])

fprintf('\n')
fprintf('Year Guliya Dunde Dasuopu QuelcS Sajama Huasc2 \n')
fprintf('%6.0f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f \n', ...
    [dec d18O]')

fprintf('\n')
fprintf('Year HCZ ACZ LCZ6 LCZ4 \n')
fprintf('%6.0f %8.4f %8.4f %8.4f %8.4f \n', [dec hcz acz lcz lcz4]')

fprintf('\n')
fprintf('Year Recon .025 .25 .5 .75 .975 \n')
fprintf('%6.0f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f \n', ... 
    [dec'; recon(1:99)' recon4(100); reconq(1:99,:)' reconq4(100,:)'])