%This simulates the reduction in error (RE) Monte Carlo Results of
%MM05 using Matlab rather than R code. 
%
%Results are sensitive to how synthetic proxy data are
%generated; three methods are given below (comment out all but the
%one you want).  Reassuringly, when the RE statistic is correctly 
%computed all cases give a 99% critical value
%less that calculated by MM05 or MBH98 for the AD1400 step (see Table
%3 in MM05).  Thus it is clear that MM05 are in error when they
%conclude that the RE value reported by MBH98 is insignificant.

clear;

%Instrumental records
Instrumental_Sparse;
inst_time=A(:,1); 
Temp_instr=A(:,2);

%North American tree rings
NOAMER;
proxy_time=proxy(:,2); 
proxy=proxy(:,3:end)';
proxynum=length(proxy(:,1));
  
pl_est_cal=find(proxy_time>=1902 & proxy_time<=1980); 
pl_est_ver=find(proxy_time<1902 & proxy_time>=1856); 
pl_obs_cal=find(inst_time>=1902 & inst_time<=1980); 
pl_obs_ver=find(inst_time<1902 & inst_time>=1856); 
  
for iter=1:1000,  %number of iterations
  for ct=1:proxynum,
    %Uncomment below lines for different methods of producing synthetic records
    %proxyrand(ct,:)=phaserandomize(proxy(ct,:)-mean(proxy(ct,:))); %tree-like auto-correlated noise
    %proxyrand(ct,:)=randn(1,length(proxy));         %white noise
    proxyrand(ct,:)=cumsum(randn(1,length(proxy))); %red noise
  end;    
  [U S x]=svds(proxyrand,1); 
  xfix=(x-mean(x(pl_est_cal)))*std(Temp_instr(pl_obs_cal))/std(x(pl_est_cal))+mean(Temp_instr(pl_obs_cal));
  xMM=(x-mean(x(pl_est_cal)))*std(Temp_instr(pl_obs_cal))/(4*std(x(pl_est_cal)))+mean(Temp_instr(pl_obs_cal));
  REfix(iter)=rePH(xfix(pl_est_ver),Temp_instr(pl_obs_ver));
  REMM(iter)=rePH(xMM(pl_est_ver),Temp_instr(pl_obs_ver));
  display([iter REMM(iter)  REfix(iter)]);
end;


REMM=sort(REMM); 
cvMM=REMM([.9 .95 .99]*length(REMM));

REfix=sort(REfix); 
cvfix=REfix([.9 .95 .99]*length(REfix));


display(['90% 95% and 99% critical value estimates for the RE statistics']);
display(['MM05  simulation: ',num2str(cvMM)]);
display(['Fixed simulation: ',num2str(cvfix)]);