clear 
randn('state',sum(100*clock)); %set seed for random number generator;
 
pp=load('c:\scozztemann3\newtemp\CRU_NH_reform'); % NHmean temp 1850-2006 CRU data
%pp=load('c:\scozztemann3\newtemp\HAD_NH_reform'); % NHmean temp 1850-2006 CRU data
%pp=load('c:\scozztemann3\newtemp\iCRU_NH_reform'); % NHmean temp 1850-2006 CRU data
%pp=load('c:\scozztemann3\newtemp\iHAD_NH_reform'); % NHmean temp 1850-2006 CRU data

%pp=load('c:\scozztemann3\newtemp\CRU_SH_reform'); % NHmean temp 1850-2006 CRU data
%pp=load('c:\scozztemann3\newtemp\HAD_SH_reform'); % NHmean temp 1850-2006 CRU data
%pp=load('c:\scozztemann3\newtemp\iCRU_SH_reform'); % NHmean temp 1850-2006 CRU data
%pp=load('c:\scozztemann3\newtemp\iHAD_SH_reform'); % NHmean temp 1850-2006 CRU data

aa=pp(1:146,2);

% here is where you set the verification target (i.e. nhactmean)
 data=aa;%respliced(:,1);%respliced(:,1);

% set the verification period
startverify=1; endverify=46;
startcalib=47; endcalib=146;

verificationseries=data;
keepverificationseries=data;

[nrows,ncols]=size(data);

% get mean and standard deviation of calibration period
means=nanmean(data(startcalib:endcalib,:));
stds=nanstd(data(startcalib:endcalib,:));
data=data-means;

% how many years of verification
numberofvalues=(nrows);

% get lag-1 autocorrelation coefficient for each data column
for i=1:ncols
[coef1(:,i),serr,xtxi,rsq(i),fval,stats,yhat,resids,covariance] = mregress(data(startcalib:endcalib-1,i),data(startcalib+1:endcalib,i),1);

% generate 'k' red noise "reconstructions"
for k=1:1000
    montetemp(1,k)=data(startcalib,i);
 
    % generage 'numberofvalues' years of red noise data
for j=2:numberofvalues+3
    montetemp(j,k)=(coef1(2,i).*montetemp(j-1,k)+randn);
    %montetemp(j,k)=(montetemp(j-1,k)-randn)./0.9;%coef1(2,i);
end
monte(1:numberofvalues,k)=montetemp(4:numberofvalues+3,k);

    %standardize the red noise and give it the mean and std of the
    %actual calibration series
    savemonte(:,k)=monte(:,k);
    montemean(k)=mean(monte(startcalib:endcalib,k));
    montestd(k)=std(monte(startcalib:endcalib,k));
    monte(:,k)=monte(:,k)-repmat(montemean(k),numberofvalues,1);
    monte(:,k)=monte(:,k)./repmat(montestd(k),numberofvalues,1);
    monte(:,k)=monte(:,k).*repmat(stds(i),numberofvalues,1);
    monte(:,k)=monte(:,k)+repmat(means(i),numberofvalues,1);

% to lowpass uncomment the following line and comment out the subsequent
% line
%verificationseries=lowpassmin(keepverificationseries,0.1);
verificationseries=keepverificationseries;
saveverificationseries=verificationseries;


% set up the red noise series for RE CE and rsq calculations

%to lowpass uncomment the following line
%monte(:,k)=lowpassmin(monte(:,k),0.1);
keepmonte(:,k)=monte(:,k);

% subtract the calibration mean from each series and calculate RE
verificationseries(:,i)=saveverificationseries(:,i)-repmat(means(i),numberofvalues,1);
remonte(:,k)=keepmonte(:,k)-repmat(means(i),numberofvalues,1);

diff=(verificationseries(startverify:endverify,i)-remonte(startverify:endverify,k));
sumsqdiff=nansum(nansum(diff.^2));
sumsq=nansum(nansum(verificationseries(startverify:endverify,i).^2)); %
re(k)=1-(sumsqdiff/sumsq);
clear  sumsq

%subtract the verification mean and calculate CE
verifmean(i)=nanmean(saveverificationseries(startverify:endverify,i));
cemonte(:,k)=keepmonte(:,k)-repmat(verifmean(i),numberofvalues,1);
ceverifyseries=saveverificationseries-repmat(verifmean(i),numberofvalues,1);
diff=(ceverifyseries(startverify:endverify,i)-cemonte(startverify:endverify,k));
%sumsqdiff=nansum(nansum(diff.^2));
sumsq=nansum(nansum(ceverifyseries(startverify:endverify,i).^2)); %
ce(k)=1-(sumsqdiff/sumsq);
clear sumsqdiff sumsq

% now get the r-squared
[coeff,S_err,XTXI,rsq(k),Fval,stats,Y_hat,resid,junk] = mregress(verificationseries(startverify:endverify,i),monte(startverify:endverify,k),1);

end
% now for each series of interest find the 50, 90,95 and 99th percentile
% for the RE Ce and Rsq
reprcentile(:,i)=prctile(re,[50 70 90 95 99])';
ceprcentile(:,i)=prctile(ce,[50 70 90 95 99])';
rsqprcentile(:,i)=prctile(rsq,[50 70 90 95 99])';

end



%###################
% for late-miss
%###################

%clear 
randn('state',sum(100*clock)); %set seed for random number generator;
 
%pp=load('c:\scozztemann3\newtemp\CRU_NH_reform'); % NHmean temp 1850-2006 CRU data
%pp=load('c:\scozztemann3\newtemp\HAD_NH_reform'); % NHmean temp 1850-2006 CRU data
%pp=load('c:\scozztemann3\newtemp\iCRU_NH_reform'); % NHmean temp 1850-2006 CRU data
%pp=load('c:\scozztemann3\newtemp\iHAD_NH_reform'); % NHmean temp 1850-2006 CRU data

%pp=load('c:\scozztemann3\newtemp\CRU_SH_reform'); % NHmean temp 1850-2006 CRU data
%pp=load('c:\scozztemann3\newtemp\HAD_SH_reform'); % NHmean temp 1850-2006 CRU data
%pp=load('c:\scozztemann3\newtemp\iCRU_SH_reform'); % NHmean temp 1850-2006 CRU data
%pp=load('c:\scozztemann3\newtemp\iHAD_SH_reform'); % NHmean temp 1850-2006 CRU data

%aa=pp(1:146,2);

% here is where you set the verification target (i.e. nhactmean)
 data=aa;%respliced(:,1);%respliced(:,1);

% set the verification period
startverify=101; endverify=146;
startcalib=1; endcalib=100;

verificationseries=data;
keepverificationseries=data;

[nrows,ncols]=size(data);

% get mean and standard deviation of calibration period
means=nanmean(data(startcalib:endcalib,:));
stds=nanstd(data(startcalib:endcalib,:));
data=data-means;

% how many years of verification
numberofvalues=(nrows);

% get lag-1 autocorrelation coefficient for each data column
for i=1:ncols
[coef1(:,i),serr,xtxi,rsq(i),fval,stats,yhat,resids,covariance] = mregress(data(startcalib:endcalib-1,i),data(startcalib+1:endcalib,i),1);

% generate 'k' red noise "reconstructions"
for k=1:1000
    montetemp(1,k)=data(startcalib,i);
 
    % generage 'numberofvalues' years of red noise data
for j=2:numberofvalues+3
    montetemp(j,k)=(coef1(2,i).*montetemp(j-1,k)+randn);
    %montetemp(j,k)=(montetemp(j-1,k)-randn)./0.9;%coef1(2,i);
end
monte(1:numberofvalues,k)=montetemp(4:numberofvalues+3,k);

    %standardize the red noise and give it the mean and std of the
    %actual calibration series
    savemonte(:,k)=monte(:,k);
    montemean(k)=mean(monte(startcalib:endcalib,k));
    montestd(k)=std(monte(startcalib:endcalib,k));
    monte(:,k)=monte(:,k)-repmat(montemean(k),numberofvalues,1);
    monte(:,k)=monte(:,k)./repmat(montestd(k),numberofvalues,1);
    monte(:,k)=monte(:,k).*repmat(stds(i),numberofvalues,1);
    monte(:,k)=monte(:,k)+repmat(means(i),numberofvalues,1);

% to lowpass uncomment the following line and comment out the subsequent
% line
%verificationseries=lowpassmin(keepverificationseries,0.1);
verificationseries=keepverificationseries;
saveverificationseries=verificationseries;


% set up the red noise series for RE CE and rsq calculations

%to lowpass uncomment the following line
%monte(:,k)=lowpassmin(monte(:,k),0.1);
keepmonte(:,k)=monte(:,k);

% subtract the calibration mean from each series and calculate RE
verificationseries(:,i)=saveverificationseries(:,i)-repmat(means(i),numberofvalues,1);
remonte(:,k)=keepmonte(:,k)-repmat(means(i),numberofvalues,1);

diff=(verificationseries(startverify:endverify,i)-remonte(startverify:endverify,k));
sumsqdiff=nansum(nansum(diff.^2));
sumsq=nansum(nansum(verificationseries(startverify:endverify,i).^2)); %
re(k)=1-(sumsqdiff/sumsq);
clear  sumsq

%subtract the verification mean and calculate CE
verifmean(i)=nanmean(saveverificationseries(startverify:endverify,i));
cemonte(:,k)=keepmonte(:,k)-repmat(verifmean(i),numberofvalues,1);
ceverifyseries=saveverificationseries-repmat(verifmean(i),numberofvalues,1);
diff=(ceverifyseries(startverify:endverify,i)-cemonte(startverify:endverify,k));
%sumsqdiff=nansum(nansum(diff.^2));
sumsq=nansum(nansum(ceverifyseries(startverify:endverify,i).^2)); %
ce(k)=1-(sumsqdiff/sumsq);
clear sumsqdiff sumsq

% now get the r-squared
[coeff,S_err,XTXI,rsq(k),Fval,stats,Y_hat,resid,junk] = mregress(verificationseries(startverify:endverify,i),monte(startverify:endverify,k),1);

end
% now for each series of interest find the 50, 90,95 and 99th percentile
% for the RE Ce and Rsq
reprcentilel(:,i)=prctile(re,[50 70 90 95 99])';
ceprcentilel(:,i)=prctile(ce,[50 70 90 95 99])';
rsqprcentilel(:,i)=prctile(rsq,[50 70 90 95 99])';

end

ddd=[99 99 99 99 99]';

re1=(reprcentile + reprcentilel)/2;
ce1=(ceprcentile + ceprcentilel)/2;
rsq1=(rsqprcentile + rsqprcentilel)/2;

res=[reprcentile reprcentilel re1 ddd ...
    ceprcentile ceprcentilel ce1 ddd ...
    rsqprcentile rsqprcentilel rsq1 ddd]




