function [sd99,sd95,sd90,sd5]=gsbTest(sL,sN,window,dR,method)

% GSBTEST  Monte Carlo Test of Moving Correlation Significance Levels
%              
%      [sd99,sd95,sd90,sd5]=gsbTest(sL,sN,window,R,method)
%
% Implementation of a Monte Carlo type simulation of moving correlation
% functions between two weakly correlated white noise series for the purposes of constructing
% significance test for moving correlation functions based on their standard deviation (Gershunov et al. 2001).
%
% Inputs:
%   sL = length of the time series used in the running correlation analysis
%   sN = number of Monte Carlo iterations to perform (should be >= 1,000)
%   window = width (in units of the time step for the data) of the correlation window
%   dR = desired overall correlation of the simulated time series (simulations will be with 1% of this value)
%   method [optional] = Which method to use to create correlated random series [default = 1], see below
%
% Correlation Methods
%   1 = Cholsky decomposition of the covariance matrix [default]
%   2 = Matrix square root
%
% Outputs:
%  sd99 = bootstrapped 99% confidence level
%  sd95 = bootstrapped 95% confidence level
%  sd90 = bootstrapped 90% confidence level
%  sd5 = bootstrapped 5% confidence level  
%
% Example:
%  > [sd99,sd95]=gsbTest(100,1000,11,0.30,2)
%
% sd99 =
%
%    0.4002
%
% sd95 =
%
%    0.3624
%           
% NOTE: Although the solution should converge for sN >= 1000, small differences (due to the correlation between
% the two artificial time series generated in the script), may still result, especially for a small number (<<10000)
% of simulations; so it is recommended that you verify convergence with several runs and be cautious of values too
% close to your chosen test statistic (particularly at 90% or 95% CI).  You can compare the above result for 'sd95'
% to Table 1 from Gershunov et al. 2001 (interseries r = 0.30, window length = 11, 95% confidence level = 0.36).
%
% References
% ----------
% Alexander Gershunov, Niklas Schneider, and Tim Barnett, 2001. Low-frequency modulation of the ENSO-Indian 
% monsoon rainfall relationship: Signal or noise?, Journal of Climate, 14: 2486-2492.
%

% Update History:
%
% November 2009 -- v1.1 now correctly points to pairwise.m
% June 2008 -- v1.0 posted online
% October 2007 -- v1.0 Initial public version
% September 2003 -- v0.1 Initial version
% 
% Function gsbTest.m written by Kevin Anchukaitis
% Subfunction build_win is a modification of the function by David Meko (dmeko@ltrr.arizona.edu)
% Subfunction pairwise is a modification of corrpair.m by David Meko (dmeko@ltrr.arizona.edu)

if (nargin<5), method = 1; end;
if (nargin<4) error('Too Few Input Parameters'); end;

% First, create a correlation matrix from the desired correlation value for the white noise series
% Then, determine the method to be used for the construction of the correlated
% white noise series
Rmat=[1 dR;dR 1];

randn('state',sum(100*clock));  % re-seed the random number generator

if method==1 %% Cholsky decomposition
    	for i=1:sN
            C=chol(Rmat);  % Cholsky decomposition of the correlation matrix
            rcheck = 1;
             while abs(rcheck) > .005*dR; % constrain within 0.5%, could be relaxed
              Z=randn(2,sL);
              X=C'*Z; X=X';
              check = corrcoef(X);
              rcheck = check(2,1) - dR;
             end 
            % now do the moving correlation function on the series in X
            x1=build_win(X(:,1),window); x2=build_win(X(:,2),window);
            s = pairwise(x1,x2);
            mcf_var(i) = var(s); mcf_std(i) = std(s);
		end
		
	elseif method==2 %% Matrix square root
		for i=1:sN
		    rcheck = 1;
             while abs(rcheck) > .005*dR; % constrain within 0.5%, could be relaxed
              X=(sqrtm(Rmat)*randn(2,sL))';
              check = corrcoef(X);
              rcheck = check(2,1) - dR;
             end
            % now do the moving correlation function on the series in X
            x1=build_win(X(:,1),window); x2=build_win(X(:,2),window);
            s = pairwise(x1,x2);
            mcf_var(i) = var(s); mcf_std(i) = std(s);
		end   

else 
error('Invalid Series Creation Method'); 
end
    
% sort the matrix of variance or standard deviation 
sorted_mcf_var = sort(mcf_var'); sorted_mcf_std = sort(mcf_std');

flag_99 = round(0.99 * sN); flag_95 = round(0.95 * sN); flag_90 = round(0.90 * sN); flag_5 = round(0.05 * sN);
sd99 = sorted_mcf_std(flag_99); sd95 = sorted_mcf_std(flag_95); sd90 = sorted_mcf_std(flag_90); sd5 =  sorted_mcf_std(flag_5);   


function X=build_win(x,m)
% Build windowed matrix from time series vector
% Adapted from code originally written by David Meko
[mx,nx]=size(x);
if nx~=1;
    error('x must be a column vector');
end;

if m>=mx;
    error('window width must be shorter than time series!');
end;

% Make index to end times
i2 = fliplr(mx:-1:m);

% Expand to matrix
A = repmat(i2,m,1);
[mA,nA]=size(A);

% increment vector
b = (0:1:(m-1))';
B=repmat(b,1,nA);

C=A-B;
C=flipud(C);
X=x(C);



function s = pairwise(X,Y)
% modified from David Meko's corrpair.m
[mX,nX]=size(X);
[mY,nY]=size(Y);
x_bar = mean(X);
X_bar = repmat(x_bar,mX,1);
y_bar = mean(Y);
Y_bar = repmat(y_bar,mY,1);

% Matrices of std devs; these are computed with "N-1" in denominator
xstd = std(X);
Xstd = repmat(xstd,mX,1);
ystd = std(Y);
Ystd = repmat(ystd,mY,1);

% Matrices of departures
Dx = X - X_bar;
Dy = Y - Y_bar;

% Matrices of squared departures
D2x = Dx .* Dx;
D2y = Dy .* Dy;

% Product of departures
D2 =  Dx .* Dy;

% Sample covariance
D = sum(D2)/(mX-1); % rv

% Std of x time std of y
xystd = xstd .* ystd; % a rv

% Correlation 
s = D ./ xystd;