%Set a flags here for hybrid.
ishybrid=1; % use 1 for hybrid 0 for nonhybrid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Load necessary data files. The setup below is all based on .mat input
%files
load csmprox00.mat; %Name of proxy data file. The variable in the file is the same as the file name.
pprox = csmprox00; %Sets the name of the proxy data file (from the previous line)
load csmsub.mat; %Same story here for the model field
csm = csmsub;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nyears,nproxies] = size(pprox);
[nyears,ngrids] = size(csm);
%Set up matrix for RegEM.
composite=zeros(nyears,ngrids+nproxies); %First initialize a matrix of appropriate size.
composite=composite./0;
composite(1:1131,1:669)= csm(1:1131,1:669); %Now put in model data.
composite(1:1131,670:669+nproxies)=pprox(1:1131,1:104);%Add the pseudo-proxies to the matrix.
composite=composite(1:1131,1:669+nproxies); %Finally, cut composite at 1980 (for 1901-1980 calibration interval).
%Standardize everything
means=nanmean(composite(1:1131,:));%This is for the model interval from 850(1) to 1980(1131)
stdevs=nanstd(composite(1:1131,:));
composite=composite-repmat(means,1131,1);%Subtract the mean.
composite=composite./repmat(stdevs,1131,1);%Divide by the standard deviations
%For M05, standardization is over total interval. Now erase pre-cal
%interval.
composite(1:1006,1:669) = NaN;
%These are necessary inputs for RegEM.
[nyears,nvars] = size(composite);
weight=[];%For weight = emptyset, RegEM uses default weighting (see line 400 in regem.m)
centerend=1006;%This is set for initial centering of the data matrix.
%This is the end of the reconstruction period. In this code the
%beginning has a default setting of 1. The data are centered from 850
%to 1856. See center.m for more.
if ishybrid==0
X=composite(1:nyears,:);
% Options for regem
OPTIONS.regress = 'mridge';
%OPTIONS.regpar = 100;
OPTIONS.stagtol = 5e-3;
OPTIONS.maxit = 200;
OPTIONS.inflation = 1;
OPTIONS.disp = 1;
OPTIONS.relvar_res = 0.050;
OPTIONS.minvarfrac=0.95;
%OPTIONS.Xcmp = keepers;
[X, M, C, Xerr, Xmis, dXmis] = regem(X, OPTIONS, weight,centerend);
elseif ishybrid==1;
cutoff=20; %years
%nproxies=112;
lowlim=1./(cutoff./2)
%start with some weights
%tricky here, a higher number is a lower weight. The actual weight in RegEM
%is the reciprocal of the weight here. i.e. RegEM devides by the weights.
clear *weight*
weight(1:nvars)=1; %thats the pcproxies
weight=diag(weight);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nrows,ncols]=size(composite)
gotstart=zeros(1,ncols);
gotstop=zeros(1,ncols);
start=zeros(1,ncols);
stop=zeros(1,ncols);
stop(:)=nrows;
[lowa,lowb]=butter(5,lowlim,'low');
%%%%%%%%%%%%%%
%%%% This section finds the starting and stopping year for each column. It assumes that
%%%% there are no missing values between start and stop, otherwise filtering chokes.
for iii=1:nrows
for jj=1:ncols
if ~isnan(composite(iii,jj)) & gotstart(jj)~=1
start(jj)=iii;
gotstart(jj)=1;
elseif isnan(composite(iii,jj)) & gotstart(jj)==1 & gotstop(jj)~=1
stop(jj)=iii-1;
gotstop(jj)=1;
end
end
end
%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Apply the zero-phase butterworth filter but only to the proxies
highf=zeros(nrows,ncols);
lowf=zeros(nrows,ncols);
lowf=lowf./lowf;
highf=highf./highf;
for k=1:ncols
clear y y2
data=composite(start(k):stop(k),k);
lowf(start(k):stop(k),k)=lowpass(data,0.05,1,1);
highf(start(k):stop(k),k)=data-lowf(start(k):stop(k),k);
end
%% Now scale lowf and highf proxies only by the percentage of variance
%each contains relative to the unfiltered series (ignore NaNs for this).
lowfpct=(nanstd(lowf).^2)./(nanstd(composite).^2);
lowfpct=1-lowfpct;
highfpct=(nanstd(highf).^2)./nanstd(composite).^2;
highfpct=1-highfpct;
highfpct(1:ncols-nproxies)=1; % keep instrumental weights=1
lowfpct(1:ncols-nproxies)=1; % keep instrumental weights 1
%The weights need to be a diagonal matrix
highweight=diag(highfpct);
highweight=highweight*weight; %combine the weights
%lowweight=diag(lowfpct);
%lowweight=lowweight*weight; %combine the weights
%%%%%%%%%%%%%%%%%%%%%%%%
%% NOW RUN THE high FREQUENCY PART
clear lowweight lowfpct lowf k j jj jjj ii iii *pct*
X=highf(1:nyears,:);
% Options for regem
OPTIONS.regress = 'mridge';
%OPTIONS.regpar = 100
OPTIONS.stagtol = 5e-2;
OPTIONS.maxit = 200;
OPTIONS.inflation = 1;
OPTIONS.disp = 1;
OPTIONS.relvar_res = 0.05;
OPTIONS.minvarfrac=0.95;
%OPTIONS.Xcmp = keepers;
%
[X, M, C, Xerr, Xmis, dXmis] = regemhigh_hpar(X,OPTIONS,highweight,centerend);
end % end the ishybrid part
%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%
%Last Section
csm00_m05_high =X(1:1131,1:669); %LHS is the name of the completed Regem reconstruction.
save finished_csm00_m05_sub_high
clear C M OPTIONS X Xerr Xmis pprox echog ngrids allcenterend centerend composite constant dXmis i ishybrid mbhstepwise means nproxies nvars nyears season stdevs temp weight yearsteps
%%%%Finished!