function [w, sig, c, c1, k] = bds (series, maxdim, distance, flag, maxram) %BDS Brock, Dechert & Scheinkman test for independence based on the correlation dimension % % [W, SIG, C, C1, K] = BDS (SERIES, MAXDIM, DISTANCE, METHOD, MAXRAM) % % uses - time-series vector SERIES (1), % - dimensional distance % * either defined as fraction DISTANCE of the standard deviation of SERIES % if FLAG = 0, % * or defined such that the one dimensional correlation integral of SERIES % is equal to DISTANCE if FLAG = 1 (2), % - not more than MAXRAM megabytes of memory for the computation (3), % % to compute - BDS statistics W for each dimension between 2 and MAXDIM (4), % - significance levels SIG at which the null hypotheses of no dependence are % rejected ASYMPTOTICALLY (use companion function BDSSIG.M for finite % samples) against (almost) any type of linear and non-linear dependence (5), % - correlation integral estimates C for each dimension M between 2 and MAXDIM, % - first-order correlation integral estimates C1 computed over the last N-M+1 % observations, and % - parameter estimate K (6). % % (1) SERIES is normally a vector of residuals obtained from a regression, but it can also % be any other stationary time series. % (2) The default settings are DISTANCE = 1.5 and FLAG = 0. The BDS statistic appears to % be most efficient estimated if the measure of dimensional distance EPSILON is chosen % such that the first-order correlation integral estimate (C1) lies around 0.7 (see % Kanzler, 1998, forthcoming). For settings DISTANCE = 0.7 and FLAG = 1, the % programme will chose EPSILON accordingly. Unfortunately, the cost of finding optimal % EPSILON is quite high in terms of CPU time and required memory. For a near-normal % distribution, the default settings achieve the same without any extra computational % burden. % (3) The default setting is MAXRAM = 150, which is recommended for a system with 192MB % physical RAM installed. The programme is highly optimised as to maximise speed given % available memory, so it is very important to specify MAXRAM correctly as the amount % of physical memory available AFTER starting MATLAB, loading any data and running % other applications concurrently. The smaller the amount of RAM available to the % programme (in relation to the length of SERIES), the slower the algorithm chosen % from six alternatives. However, if MAXRAM is chosen too large, MATLAB will make use % of virtual (hard-disk) memory, and this will slow down computation considerably. % (4) The default setting for MAXDIM is 2. For MAXDIM = 1, W and SIG are empty. % (5) A vector of NaN is returned if the MATLAB Statistics Toolbox is not installed. % (6) The BDS statistic W(M) is a function of C(M), C1(1), C1(M) and K, and these % estimates are normally of no further interest. % % See Kanzler (1998) for some explanation of the main parts of the algorithm (other % explanations are commented into the below code), for a detailed investigation of the % finite-sample properties of the BDS statistic, for tables of small-sample quantiles and % for a comparison with software by Dechert (1988) and LeBaron (1988, 1990, 1997a, 1997b). % These and other important references can be found at the end of the script. % % * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * % * All rights reserved. This script may be redistributed if it is left unaltered in * % * its entirety (619 lines, 31422 bytes) and if nothing is charged for redistribution. * % * Usage of the programme in applications and alterations of the code should be * % * referenced properly. See http://users.ox.ac.uk/~econlrk for updated versions. * % * The author appreciates suggestions for improvement or other feedback. * % * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * % % Copyright (c) 14 April 1998 by Ludwig Kanzler % Department of Economics, University of Oxford % Postal: Christ Church, Oxford OX1 1DP, England % E-mail: ludwig.kanzler@economics.oxford.ac.uk % $ Revision: 2.41 $ $ Date: 15 September 1998 $ % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % Executable part of main function BDS.M starts here % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %%%%%%%%%%%%%%%%%%%%%% Check and transformation of input arguments %%%%%%%%%%%%%%%%%%%%%% if nargin < 5 maxram = 150; elseif maxram > 500 disp('Are you sure you have so much memory available?') error('If so, you need to edit the code, otherwise try again with a lower value.') end if nargin < 4 flag = 0; elseif ~any(flag == [0 1]) error('Unknown method for determining dimensional distance; try again with 0 or 1.') end if nargin < 3 distance = 1.5; elseif distance < 0 error('The dimensional distance parameter must be positive.') elseif flag == 1 & distance > 1 error('The correlation integral cannot exceed 1.') end if nargin < 2 maxdim = 2; elseif maxdim < 1 error('The dimension needs to be at least 1.'); end if nargin < 1 error('Cannot compute the BDS statistic on nothing.') end [rows,cols] = size(series); if rows > 1 & cols == 1 n = rows; series = series'; elseif cols > 1 & rows == 1 n = cols; elseif cols > 1 & rows > 1 n = cols*rows; series = series(:)'; % transformation into a row vector disp(sprintf('\aTransformed matrix input into a single column.')) else error('Cannot compute the BDS statistic on a scalar!') end %%%%%%%%%%%% Determination of and preparations for fastest method given MAXRAM %%%%%%%%%%% fastbuild = 0.000016 * (1:52) .* pow2(1:52); % memory requirements slowbuild = 0.000045 * pow2(1:52); % for the various holdinfo = 0.000005 * pow2(1:52); % algorithms in wordtable = 0.000008 * n^2 ./ (1:52); % megabytes for bitandop = 0.000024 * n^2 ./ (1:52); % given N [ram1, bits1] = min(fastbuild + holdinfo + wordtable + bitandop); % number of bits for [ram2, bits2] = min(fastbuild + holdinfo + wordtable); % which each of six [ram3, bits3] = min(slowbuild + holdinfo + wordtable + bitandop); % methods uses minimum [ram4, bits4] = min(slowbuild + holdinfo + wordtable); % memory; this memory [ram5, bits5] = min( wordtable + bitandop); % is given by [ram6, bits6] = min( wordtable); % ram1, ram2,..., ram6 if ram1 < maxram | ram2 < maxram if ram1 < maxram method = 1; bits = bits1; ram = ram1; else method = 2; % maximum number of rows to put bits = bits2; ram = ram2; % through BITAND and bit-counting stepping = floor((maxram-ram)*bits/n/0.000024); % algorithm without exceeding MAXRAM end % Vector BITINFO lists the number of bits set for each integer between 0 and 2^bits % (corresponding to the indices of the vector shifted by 1). See Kanzler (1998) for an % explanation. bitinfo = uint8(sum(rem(floor((0:pow2(bits)-1)'*pow2(1-bits:0)),2),2)); elseif ram3 < maxram | ram4 < maxram if ram3 < maxram method = 3; bits = bits3; ram = ram3; else method = 4; bits = bits4; ram = ram4; stepping = floor((maxram - ram) * bits / n / 0.000024); end bitinfo(1:pow2(bits), :) = uint8(0); % the same as above, but created through for bit = 1 : bits % a loop, which consumes less memory bitinfo(1:pow2(bits)) = sum([bitinfo, ... kron(ones(pow2(bits-bit),1), [zeros(pow2(bit-1),1); ones(pow2(bit-1),1)])],2); end elseif ram5 < maxram | ram6 < maxram if ram5 < maxram method = 5; bits = bits5; ram = ram5; else method = 6; bits = bits6; ram = ram6; stepping = floor((maxram - ram) * bits / n / 0.000024); end else disp('Insufficient amount of memory. Allocate more memory to the system') disp('or reduce the number of observations, then try again.') error(' ') end %%%%%%%%%%%%%%%%%%%%% Determination of dimensional distance EPSILON %%%%%%%%%%%%%%%%%%%%%% % The empirical investigation by Kanzler (1998) shows that choosing EPSILON such that the % first-order correlation integral is around 0.7 yields the most efficient estimation of % low-dimensional BDS statistics. Hence the objective here is to choose EPSILON such that, % say, 70% of all observations lie within distance EPSILON to each other. If desired, the % programme first determines EPSILON as to fulfil this or a similar requirement. % % The conceptually simplest way of setting up the calculation of distance among all % observations is to define a two-dimensional table D (for "distance") of length and width % N and assign to each co-ordinate (x,y) the result of the problem ABS(x-y). % % In principle, the entire table could thus be created with the following one-line % statement: % D = ABS( SERIES(ONES(1000,1),:)' - SERIES(ONES(1000,1),:) ) % % Since the lower triangle of the table only replicates the upper triangle and since the % diagonal values represent own values (ones) which are not desired to be included in the % calculation, only the upper triangle receives further attention. % % Unfortunately, sewing all the row vectors of the upper triangle together to form one % single (row) vector makes indexing very messy. To aid understanding of the vector-space % indexing used here (as well as in the optional sub-function further below), one may wish % to refer to the following exemplary matrix table (N=7): % % Using this example, it is easy to verify % * * * *c o l u m n* * * * that column vector I is defined by the % I 1 2 3 4 5 6 7 following indices in vector space: % I+(0 : I-2)*N - CUMSUM(1 : I-1) % * 1 * 1 2 3 4 5 6 % * 2 . * 7 8 9 10 11 More generally, column vector I starting only % r 3 . . * 12 13 14 15 in row J is: % o 4 . . . * 16 17 18 I+(J-1 : I-2)*N - SUM(1:J-1)-CUMSUM(J : I-1) % w 5 . . . . * 19 20 % * 6 . . . . . * 21 Row vector I is given by indices: % * 7 . . . . . . * 1+(I-1)*(N-1)-SUM(1:I-2) : I*(N-1)-SUM(1:I-1) % % (A formal derivation of the above formulae is beyond the scope of this script.) % % To calculate a percentile of the distribution of distance values, the row vector is % sorted (unfortunately, this requires a lot of time and RAM in MATLAB). if ~flag demeaned = series-sum(series)/n; % fastest algorithm for epsilon = distance * sqrt(demeaned*demeaned'/(n-1)); % computing the standard clear demeaned % to save memory % deviation of SERIES elseif 0.000008 * 3 * sum(1:n-1) < maxram % check memory requirements for DIST and sorting dist(1:sum(1:n-1)) = 0; for i = 1 : n-1 dist(1+(i-1)*(n-1)-sum(0:i-2):i*(n-1)-sum(1:i-1)) = abs(series(i+1:n)-series(i)); end sorted = sort(dist); epsilon = sorted(round(distance*sum(1:n-1))); % DISTANCEth percentile of SORTED series clear dist sorted else error('Insufficient RAM to compute EPSILON; allocate more memory or use METHOD = 1.') end %%%%%%%%%%%% Computation and storage of one-dimensional distance information %%%%%%%%%%%% % Similarly to the above, a two-dimensional table C (for "close") of length and width N % can be defined by assigning to each co-ordinate (x,y) the result of the problem % ABS(x-y) <= EPSILON; (x,y) assumes the value 1 if the statement is true and 0 otherwise. % % Formally, for given EPSILON: % C(x,y) = 1 if ABS(x-y) <= EPSILON % = 0 otherwise % % Once again, the resulting information needs to be stored in the most efficient way. % In this implementation, this is done by chopping each row of the table into "words" of % several bits, the precise number of bits per word being determined by the above % algorithms. One "word" is thus represented by one integer. This slashes the size of the % table by the number of bits. See Kanzler (1998) for more details. % % The below routine stores all rows of the upper triangle of the conceptual table % (described in Kanzler, 1998) left-aligned and assigns zeros to all other elements. % % As will also be explained further below, the computation of parameter K requires the sum % of each FULL row, i.e. each row including the elements in the lower triangle and on the % diagonal. The "missing" bits correspond to the sums over each column in the upper % triangle, and these sums are also computed and stored in the below loop. And to make % matters simple, diagonal values are allocated to the column sums by initialising them % with value 1. See also Kanzler (1998). colsum(1:n) = 1; rowsum(1:n) = 0; nwords = ceil((n-1)/bits); wrdmtrx(1:n-1,1:nwords) = 0; % initialisation of bit-word table for row = 1 : n-1 bitvec = abs(series(1+row:n) - series(row)) <= epsilon; rowsum(row) = sum(bitvec); colsum(1+row:n) = colsum(1+row:n) + bitvec; nwords = ceil((n-row)/bits); wrdmtrx(row,1:nwords) = (reshape([bitvec,zeros(1,nwords*bits-n+row)],...%transformation bits, nwords)' *pow2(0:bits-1)')'; %into bit-words end clear series bitvec %%%%%%%%%%%%%%%%%% Computation of one-dimensional correlation estimates %%%%%%%%%%%%%%%%%% % C1(1), the fraction (or estimated probability) of pairs in SERIES being "close" in the % first dimension is just the average over ALL unique elements. C1(1) is hence the most % efficient estimator of C(1), and the resulting estimate is used in the computation of % SIGMA(M) further below. % % However, for the difference term C(M) - C(1)^M of the BDS statistic (see further below) % to follow SIGMA asymptotically, both C(M) and C(1) need to be estimated over the same % length vector, and so MAXDIM different C1's need to be estimated here: % % N N % C1(M) = 2/(N-M+1)/(N-M) * SUM SUM B(S,T) % S=M T=S+1 % % Each C1(M) is easily computed from the sum of all bits set in rows M to N-1 divided by % the appropriate total number of bits. bitsum(maxdim:-1:1) = cumsum([sum(rowsum(maxdim:n-1)), rowsum(maxdim-1:-1:1)]); c1 (maxdim:-1:1) = bitsum(maxdim:-1:1) ./ cumsum([sum(1:n-maxdim), n-maxdim+1 : n-1]); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Computation of parameter K %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % A parameter needed to estimate SIGMA(M) is K, which is defined as: % % N N N % K = 6/N/(N-1)/(N-2)* SUM SUM SUM {C(T,S)*C(S,R) + C(T,R)*C(R,S) + C(S,T)*C(T,R)} / 3 % T=1 T=T+1 R=S+1 % % As is readily apparent, a literary computation of the above would be very processing % intensive, e.g.: % HT(1 : N) = 0; % FOR T = 1 : N % HS(1 : N-T) = 0; % FOR S = T+1 : N % HR(1 : N-S) = 0; % FOR R = S + 1 : N % HR(R-S) = (C(T,S)*C(S,R) + C(T,R)*C(R,S) + C(S,T)*C(T,R)) / 3; % END % HS(S-T) = SUM(HR(1 : N-S)); % END % HT(T)= SUM(HS(1 : N-T)); % END % K = SUM(HT) * 6 / N/(N-1)/(N-2); % % To understand what k actually estimates, and how this estimation can be made % computationally more efficient, see Kanzler (1998). % % The above FOR loop computes the sum over each row and over each column including the % diagonal in the upper triangle. To compute K from this, the sum of the squares of the % row and column sums needs to be adjusted as reasoned above, whereby the sum of all % elements in table C is given by twice the sum of all vector elements plus the diagonal % values. fullsum = rowsum + colsum; k = (fullsum*fullsum' + 2*n - 3*(2*bitsum(1)+n)) / n/(n-1)/(n-2); clear rowsum colsum fullsum bitsum %%%%%%%%%% Computation of correlation estimates and SIGMA for higher dimensions %%%%%%%%%% % C(M), the M-dimension correlation % estimate, is defined as: N N M-1 % C(M) = 2/(N-M+1)/(N-M) * SUM SUM PROD B(S-J, T-J) % S=M T=S+1 J=0 % % To see how C can be computed for M > 1, see Kanzler (1998). % % In practice, the required BITAND-operation can be performed on the entire table at once % by replacing the entire table between rows M and N-1 with the result of the BITAND- % operation between the table formed by rows M to N-1 and M-1 to N-2. But this works only % if sufficient memory is available (methods 1, 3 and 5). Otherwise, the BITAND-operation % has to be performed by looping BACKWARDS through the table, taking as many rows as % possible at once (methods 2, 4 and 6). % % The number of bits set in rows M to N-1 (inclusive) is counted either in one go or % through the above loop by looking up the number of bits set for each integer (LeBaron, % 1997, uses a similar method), or, if memory was insufficient to create the required % BITINFO array, by column-wise brute-force counting. % % MATLAB uses logarithms to compute powers, and this can result in minute deficiencies in % accuracy. To avoid this, integer powers are computed by separate functions in this % script (see further below). Otherwise, SIGMA would be calculated as follows: % sigma(m-1) = 2*sqrt(k^m + 2*k.^(m-(1:m-1))*(c1(1).^(2*(1:m-1)))'... % + (m-1)^2*c1(1)^(2*m) - m^2*k*c1(1)^(2*m-2)); for m = 2 : maxdim bitcount = 0; if sum(method == [1 3]) wrdmtrx(m:n-1,:) = bitand(wrdmtrx(m:n-1,:),wrdmtrx(m-1:n-2,:)); % BITAND and bit bitcount = sum(sum(bitinfo(wrdmtrx(m:n-1,:)+1))); % count all at once elseif sum(method == [2 4]) for row = n-stepping : -stepping : m+1 % BITAND wrdmtrx(row:row+stepping-1,:) = bitand(wrdmtrx(row:... % and bit row+stepping-1,:), wrdmtrx(row-1:row+stepping-2,:)); % count in bitcount=bitcount+sum(sum(bitinfo(wrdmtrx(row:row+stepping-1,:)+1))); % backward end % loops wrdmtrx(m:row-1,:) = bitand(wrdmtrx(m:row-1,:), wrdmtrx(m-1:row-2,:)); % through bitcount = bitcount + sum(sum(bitinfo(wrdmtrx(m:row-1,:)+1))); % the table elseif method == 5 wrdmtrx(m:n-1,:) = bitand(wrdmtrx(m:n-1,:), wrdmtrx(m-1:n-2,:)); % BITAND at once... for col = 1 : ceil((n-1)/bits) % bit count bitcount = bitcount + sum(sum(rem(floor(wrdmtrx(m:... % by brute force n-1-(col-1)*bits, col) * pow2(1-bits:0)), 2))); % in loops end else for row = n-stepping : -stepping : m+1 wrdmtrx(row:row+stepping-1,:) = bitand(wrdmtrx(row:... % BITAND row+stepping-1,:), wrdmtrx(row-1:row+stepping-2,:)); % opera- end % tions wrdmtrx(m:row-1,:) = bitand(wrdmtrx(m:row-1,:), wrdmtrx(m-1:row-2,:)); % and brute- for col = 1 : ceil((n-1)/bits) % force bit bitcount = bitcount + sum(sum(rem(floor(wrdmtrx(m:... % counting n-1-(col-1)*bits, col) * pow2(1-bits:0)),2))); % in loops end end c(m-1) = bitcount / sum(1:n-m); % indexing of sigma(m-1) = 2*sqrt(prod(ones(1,m)*k) + 2*ivp(k,m-(1:m-1),m-1)... % C and SIGMA *(ivp(c1(1),2*(1:m-1),m-1))' + (m-1)*(m-1)... % runs from 1 *prod(ones(1,2*m)*c1(1)) - m*m*k*prod(ones(1,2*m-2)*c1(1))); % to MAXDIM-1 end clear wrdmtrx %%%%%%%%%%%%%%% Computation of the BDS statistic and level of significance %%%%%%%%%%%%%%% % Under the null hypothesis of independence, it is obvious that the time-series process % has the property C(1)^M = C(M). In finite samples, C(1) and C(M) are consistently % estimated by C1(M) and C(M) as above. Also, Brock et al. (1996) show that the standard % deviation of the difference C(M) - C1(M)^M can be consistently estimated by SIGMA(M) % divided by SQRT(N-M+1), where: % % M-1 % SIGMA(M)^2 = 4* [K^M + 2* SUM {K^(M-J)* C^(2*J)} + (M-1)^2* C^(2*M) - M^2* K* C^(2*M-2)] % J=1 % % and C = C1(1) and K as above. % % For given N and EPSILON, the BDS Statistic C(M) - C1(M)^M % is defined as the ratio of the two terms: W(M) = SQRT(N-M+1) * -------------- % SIGMA(M) % % Since it follows asymptotically the normal distribution with mean 0 and variance 1, % hypothesis testing is straightforward. If available, this is done here using function % NORMCDF of the MATLAB Statistics Toolbox. % % Integer powers are again calculated by a sub-routine which is more accurate than the % MATLAB built-in power function; without using the sub-routine, the line for calculating % W would be: w = sqrt(n-(2:maxdim)+1) .* (c - c1(2:maxdim).^(2:maxdim)) ./ sigma; if maxdim > 1 w = sqrt(n-(2:maxdim)+1) .* (c - idvp(c1(2:maxdim), 2:maxdim, maxdim-1)) ./ sigma; if exist('normcdf.m','file') & nargout > 1 sig = min(normcdf(w,0,1), 1-normcdf(w,0,1)) * 2; elseif nargout > 1 sig(1:maxdim-1) = NaN; end else w = []; sig = []; c = []; end %%%%%%%%%%%%%%%%%%%%%%% Sub-functions for computing integer powers %%%%%%%%%%%%%%%%%%%%%%% function ipow = ivp (base, intpowvec, veclen) ipow(1 : veclen) = 0; for j = 1 : veclen ipow(j) = prod(ones(1, intpowvec(j)) * base); end function ipow = idvp (basevec, intpowvec, veclen) ipow(1 : veclen) = 0; for j = 1 : veclen ipow(j) = prod(ones(1, intpowvec(j)) * basevec(j)); end % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % Executable part of main function BDS.M ends here % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % The following sub-function is not actually used by the main function and only % % included for the benefit of those who would like to implement the BDS test in a % % language which is either incapable of or inefficient in handling bit-wise AND- % % operations, or those who would like to cross-check the above computation. Deleting % % the sub-function from the script will NOT result in any increase in performance. % % % % To use the function, save the remainder of this code in a file named BDSNOBIT.M. % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % function w = bdsnobit (series, maxdim, eps) %BDSNOBIT BDS test for independence IMPLEMENTED WITHOUT USING BIT-WISE FUNCTIONS % % Only such comments which relate exclusively to this implementation of the test and % which cannot be found in the main function are included below. % % Copyright (c) 14 April 1998 by Ludwig Kanzler % Department of Economics, University of Oxford % Postal: Christ Church, Oxford OX1 1DP, England % E-mail: ludwig.kanzler@economics.oxford.ac.uk % $ Revision: 1.3 $ $ Date: 30 April 1998 $ %%%%%%%%%%%%%%%%%%%%%% Check and transformation of input arguments %%%%%%%%%%%%%%%%%%%%%% if nargin < 3 eps = 1; if nargin == 1 maxdim = 2; elseif maxdim < 2 error('MAXDIM needs to be at least 2!'); end end epsilon = std(series)*eps; series = series(:)'; % PAIRS is the total number of unique pairs which can be n = length(series); % formed from all observations (note that while this is pairs = sum(1:n-1); % just (N-1)*N/2, MATLAB computes SUM(1:N-1) twice as fast!) %%%%%%%%%%%% Computation and storage of one-dimensional distance information %%%%%%%%%%%% % Recall that in the implementation of the main function above, table C is stored in bit- % representation. When this is not possible or desirable, the second best method is to use % one continuous vector of unassigned 8-bit integers (called UINT8). This, however, % requires version 5.1 or higher, and a similar option may not be available in other high- % level languages. Implementation does not depend on the ability to use unassigned low-bit % integers and would work equally with double-precision integers, but the memory % requirements would, of course, be higher. Using UINT8's is still a rather inefficient % way of storing zeros and ones, which in principle require only a single bit each. On the % PC, MATLAB actually requires "only" around 5 bytes for each UNIT8. b(1:pairs) = uint8(0); for i = 1 : n-1 b(1+(i-1)*(n-1)-sum(0:i-2):i*(n-1)-sum(1:i-1)) = abs(series(i+1:n)-series(i))<=epsilon; end clear series %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Computation of parameter K %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% sums(1 : n) = 0; for i = 1 : n sums(i) = sum(b(i+(0 : i-2)*n - cumsum(1 : i-1)))... % sum over column I + 1 ... % diagonal element + sum(b(1+(i-1)*(n-1)-sum(1:i-2) : i*(n-1)-sum(1:i-1))); % sum over row I end k = (sum(sums.^2) + 2*n - 3*(2*sum(b)+n)) / n/(n-1)/(n-2); %%%%%%%%%%%%%%%%%% Computation of one-dimensional correlation estimates %%%%%%%%%%%%%%%%%% bitsum(1:maxdim) = sum(b(1+(maxdim-1)*(n-1)-sum(0:maxdim-2) : pairs)); for m = maxdim-1 : -1 : 1 bitsum(m) = bitsum(m+1) + sum(b(1+(m-1)*(n-1)-sum(0:m-2):m*(n-1)-sum(1:m-1))); end c1(maxdim:-1:1) = bitsum(maxdim:-1:1) ./ cumsum([sum(1:n-maxdim), n-maxdim+1 : n-1]); %%%%%%%%%% Computation of correlation estimates and SIGMA for higher dimensions %%%%%%%%%% for m = 2 : maxdim % Indexing in vector space once again follows the rules set out above. Multiplication % is done by moving up column by column into north-west direction, so counter I runs % backwards in the below WHILE loop until the Mth column (from the left) is reached: i = n; while i - m % Multiplication is not defined on UINT8 variables and translating the columns % twice, once from UINT8 to DOUBLE integer and then back to UINT8, would be % inefficient, so it is better to sum entries (this operation - undocumented by % MATLAB - is defined, and even faster than the documented FIND function!) and % compare them against the value 2: b(i + (m-1 : i-2)*n - sum(1:m-1) - cumsum(m : i-1)) = ... sum([ b(i + (m-1 : i-2)*n - sum(1:m-1) - cumsum(m : i-1)); ... b(i-1 + (m-2 : i-3)*n - sum(1:m-2) - cumsum(m-1 : i-2)) ]) == 2; % The sum over each column is computed immediately after that column has been % updated. To store the column sums, the vector SUMS already used above for the row % sums is recycled (this is more memory-efficient than clearing the above SUMS % vector and defining a new vector of the column sums, because in the latter case, % MATLAB's memory space will end up being fragmented by variables K and C added to % the memory in the meantime!): sums(i) = sum(b(i + (m-1 : i-2)*n - sum(1:m-1) - cumsum(m : i-1))); i = i - 1; end c(m-1) = sum(sums(m+1:n)) / sum(1:n-m); sigma(m-1) = 2*sqrt(k^m + 2*k.^(m-(1:m-1))*(c1(1).^(2*(1:m-1)))'... % could use above + (m-1)^2*c1(1)^(2*m) - m^2*k*c1(1)^(2*m-2)); % inter-power sub- end % functions instead %%%%%%%%%%%%%%% Computation of the BDS statistic and level of significance %%%%%%%%%%%%%%% w = sqrt(n-(2:maxdim)+1) .* (c-c1(2:maxdim).^(2:maxdim)) ./ sigma; % or use sub-functions % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % Sub-function BDSNOBIT.M ends here % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % REFERENCES: % % Brock, William, Davis Dechert & José Scheinkman (1987), "A Test for Independence % Based on the Correlation Dimension", University of Wisconsin-Madison, Social % Science Research Working Paper, no. 8762 % % Brock, William, Davis Dechert, José Scheinkman & Blake LeBaron (1996), "A test for % independence based on the correlation dimension", Econometric Reviews, vol. 15, % no. 3 (August), pp. 197-235, revised version of Brock et al. (1987) % % Dechert, Davis (1988), "BDS STATS: A Program to Calculate the Statistics of the % Grassberger-Procaccia Correlation Dimension Based on the Paper "A Test for % Independence" by W. A. Brock, W. D. Dechert and J. A. Scheinkman", version 8.21 % (latest), MS-DOS software available on gopher.econ.wisc.edu % % Kanzler, Ludwig (1998), "Very Fast and Correctly Sized Estimation of the BDS Statistic", % Oxford University, Department of Economics, working paper, available on % http://users.ox.ac.uk/~econlrk % % LeBaron, Blake (1988, 1990, 1997a), "BDSTEST.C", version June 1997 (latest), C source % code available on gopher.econ.wisc.edu % % LeBaron, Blake (1997b), "A Fast Algorithm for the BDS Statistic", Studies in % Nonlinear Dynamics and Econometrics, vol. 2, pp. 53-59 % ACKNOWLEDGEMENT: % % I am grateful to Blake LeBaron for giving me the exclusive opportunity to beta-test his % C programme in its compiled version for MATLAB 5 and thus enabling me to compare the two % programmes directly. I have benefited from the many associated discussions. % End of file.