bds1.m

Neural Network in Finance (神经网络在金融界:赢得预言性的优势)全部原码。内容包括预测与估计

% 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:  [e-mail]
%                    $ 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