%% make 100 AR(1) time series

N=250;
T=200;

dats=NaN(N, T);

alp=0.6;
sig=1;



%% Generate Figures 1 and 2. 

%252 time steps:
T=252;

% there are 839 time series in the CRU annual mean temperature data set,
% but the number of observations corresponds to only about 459 complete 
% time series. (Note that this number doesn't matter for the 
% analytic calculation).  


N=459;


dats=NaN(N, T);

%Use the parameters for the BARCAST analysis of the anomalies with respect
%to 1850--2009: we want the parameters for the long reference interval
alp=0.3430;
sig=sqrt(0.4698);




% calculate the covariance matrix:
Omega=NaN(T);
for ii=1:1:T
    for jj=1:1:T
      Omega(ii,jj)=sig^2*alp^(abs(ii-jj))/(1-alp^2);
    end
end

      

dats(:,1)=randn(N,1)*sqrt(sig^2/(1-alp^2));

for kk=2:1:T
    dats(:,kk)=alp*dats(:,kk-1)+sig*randn(N,1);
end

%% plot
fill_trip=[0.7, 0.7, 0.95];
fig=figure(5); clf
set(fig, 'Position' , [10 10 1000 650], 'renderer', 'painters')
s1=subplot(2,1,1)
set(gca, 'fontsize', 14)
for kk=1:1:N
    plot(1:1:T, dats(kk,:), 'k-', 'linewidth', 1.5, 'color', [1,1,1]*0.7), hold on
end
spc=4;
boxplot(dats(:, 1:spc:end), 'colors', 'k', 'symbol', '+', 'positions', 1:spc:T)
set(gca, 'xtick', 0:20:T,  'xticklabelmode', 'auto')

plot(1:1:T, zeros(T,1), 'k-', 'linewidth', 1.5, 'color', [1,1,1]*1.0), hold on
plot(1:1:T, mean(dats), 'r-', 'linewidth', 2.5), hold on
axis([0 T+1 -3.25 3.25])
ylabel('Time Series')




std_dats=std(dats);
subplot(2,1,2)
set(gca, 'fontsize', 14)
plot(1:1:T, sqrt(Omega(1,1))*ones(1,T), 'linewidth', 1.5, 'color', [1,1,1]*0.5), hold on
plot(1:1:T, std_dats, 'k-', 'linewidth', 2.5), hold on
axis([0 T+1 0.65 0.825])
ylabel('Standard Dev.')

%% now remove from means from a 30 time step central interval:


length_int=30
Midder=mean(1:1:253);
inds_lims=floor(Midder+[-1, 1]*(length_int-1)/2);

dats_means_lims=mean(dats(:, inds_lims(1):1:inds_lims(2)), 2);
dats_adj=dats-kron(dats_means_lims, ones(1, T));


%% get the covariance matrix of the anomalies:
% A matrix:
AA=zeros(T);
AA(:, inds_lims(1):1:inds_lims(2))=-1/length(inds_lims(1):1:inds_lims(2));

Omega_Anom = (eye(T)+AA)*Omega*(eye(T)+AA)';

%%
fig=figure(6); clf
set(fig, 'Position' , [10 10 1000 650], 'renderer', 'painters')
subplot(2,1,1)
set(gca, 'fontsize', 14)
fill([inds_lims(1) inds_lims(2) inds_lims(2) inds_lims(1)], [-6 -6 6 6], fill_trip, 'linewidth', 2, 'edgecolor', fill_trip), hold on

for kk=1:1:N
    plot(1:1:T, dats_adj(kk,:), 'k-', 'linewidth', 1.5, 'color', [1,1,1]*0.7), hold on
end
boxplot(dats_adj(:, 1:spc:end), 'colors', 'k', 'symbol', '+', 'positions', 1:spc:T)
set(gca, 'xtick', 0:20:T,  'xticklabelmode', 'auto')

plot(1:1:T, zeros(T,1), 'k-', 'linewidth', 1.5, 'color', [1,1,1]*1.0), hold on
plot(1:1:T, mean(dats_adj), 'r-', 'linewidth', 2.5), hold on
axis([0 T+1 -3.25 3.25])
ylabel('Time Series')

std_dats_adj=std(dats_adj);
subplot(2,1,2)
fill([inds_lims(1) inds_lims(2) inds_lims(2) inds_lims(1)], [-6 -6 6 6], fill_trip, 'linewidth', 2, 'edgecolor', fill_trip), hold on
set(gca, 'fontsize', 14)
plot(1:1:T, sqrt(Omega(1,1))*ones(1,T), 'linewidth',  1.5, 'color', [1,1,1]*0.5), hold on
plot(1:1:T, std_dats_adj, 'k-', 'linewidth', 2.5), hold on
plot(1:1:T, diag(Omega_Anom).^(1/2), 'k-', 'linewidth',  1.5, 'color', [0,0,1]*1.0), hold on
axis([0 T+1 0.65 0.825])
ylabel('Standard Dev.')



