qc2 = qc;
for jj=1:length(qc);
temp1 = find( (X<=quantile(X,1-qc(jj))).*(Y<=quantile(Y,1-qc(jj))) );
temp2 = find( (X>=quantile(X,qc(jj))).*(Y>=quantile(Y,qc(jj))) );
if ((length(temp1)<3)+(length(temp2)<3))>0 % then not enough data to compute a correlation in at least one tail
qc2 = setdiff(qc2,qc(jj));
end
end
% now I am sure that each cut-off has enough data
qc2 = sortrows(unique([qc2;1-qc2])); % sorting and making sure that 0 does not appear twice
if isempty(find(qc==0.5))==0 % then q=0.5 is a cut-off
qc1 = [qc2(1:(end+1)/2);qc2((end+1)/2:end)];
else
qc1 = qc2;
end
q = qc1(end/2+1:end);
k = length(q);
xi = nines(T,k); % variable used in constructing test statistic
rhodiff = nines(k,1);
for jj = 1:k;
temp1 = find( (X<=quantile(X,1-qc(jj))).*(Y<=quantile(Y,1-qc(jj))) );
temp2 = find( (X>=quantile(X,qc(jj))).*(Y>=quantile(Y,qc(jj))) );
out1(jj,1) = corrcoef12(X(temp1),Y(temp1));
out1(end+1-jj,1) = corrcoef12(X(temp2),Y(temp2));
rhodiff(jj) = corrcoef12(X(temp2),Y(temp2)) - corrcoef12(X(temp1),Y(temp1)); % computing the test stat
x1plus = (X-mean(X(temp2)))./std(X(temp2));
x2plus = (Y-mean(Y(temp2)))./std(Y(temp2));
x1minus = (X-mean(X(temp1)))./std(X(temp1));
x2minus = (Y-mean(Y(temp1)))./std(Y(temp1));
xi(:,jj) = T/size(temp2,1)*(x1plus.*x2plus - corrcoef12(X(temp2),Y(temp2))).*( (X>=quantile(X,qc(jj))).*(Y>=quantile(Y,qc(jj))) ) - ...
T/size(temp1,1)*(x1minus.*x2minus - corrcoef12(X(temp1),Y(temp1))).*( (X<=quantile(X,1-qc(jj))).*(Y<=quantile(Y,1-qc(jj))) );
end
omegahat = newey_west(xi); % computing the covariance matrix
teststat = T*rhodiff'*inv(omegahat)*rhodiff;
pval = 1-chi2cdf(teststat,k);


雷达卡



京公网安备 11010802022788号







