These are programs for my Journal of Econometrics paper "On the distribution 
of the sample autocorrelation coefficients" (with Xiaolu Wang)
If you have questions, comments, or bug reports, please send them to 
kan@chass.utoronto.ca

Raymond Kan
Rotman School of Management
University of Toronto

Note 1: The Fortran subdirectory contains the source codes for a Fortran 
program NDLAED9.F that is called by rankone.m.  You can run mex1.m to produce 
the mex file.  The Fortran codes are based on the LAPACK 3.0 codes contributed
by Ren-Cang Li.  See the following paper for details:

"Solving secular equations stably and efficiently," Technical Report, 
Department of Mathematics, University of California, Berkeley, CA, USA. LAPACK 
working note 89. (available at http://www.netlib.org/lapack/lawnspdf/lawn89.pdf)

For cross-platform flexibility (i.e., for both 32 bit and 64 bit Windows,
I use Matlab's preprocessor macro by calling a Fortran header file using the line 

#include "fintrf.h"

at the beginning.  If your Fortran compiler does not support preprocessor macro,
then remove this line and replace mwpointer and mwsize in the declaration by
INTEGER*4 (for 32 bit Windows) or INTEGER*8 (for 64 bit Windows).

Note 2: If you receive an error message when running some of the programs,  
you may need to first install Microsoft Visual C++ 2008 SP1 Redistributable
Package (x86) (or x64 if you use 64 bit Windows).  For the 32 bit one, it is 
available at

http://www.microsoft.com/downloads/details.aspx?familyid=A5C84275-3B97-4AB7-A40D-3802B2AF5FC2&displaylang=en

For the 64 bit one, it is available at

http://www.microsoft.com/downloads/details.aspx?familyid=BA9257CA-337F-4B40-8C14-157CFDFFEE4E&displaylang=en


Version 1.0: 4/28/2009, initial release
Version 1.1: 11/14/2009, change the Fortran programs to work with
             both 32 bit and 64 bit Windows.  
Version 1.2: 11/21/2011, add the program to compute the density function
             of \hat\rho(k) and linear combination of \hat\rho(k) based on 
             the algorithm of Broda and Paolella (2009).
             Update the cdf and pdf programs to work with matrix input.

The following programs assume the length of times series is n, and the
data are i.i.d. normally distributed with the same mean and variance.

Eigenvalues of B_k:
eigpapb.m: A Matlab program that computes the eigenvalues of B_k.
           Usage: [y,r] = eigpapb(n,k);
           y is the unique eigenvalues of B_k, and r is the corresponding 
           multiplicities.
           The program calls rankone.m (which in turns calls rank1up.mex32
           or rank1up.mex64, you will need to mex your own if you are not 
           using 32 bit or 64 bit Windows system, the Fortran source code 
           is under the Fortran subdirectory).
eigslow.m: A slower version of eigpapb.m that relies on Matlab's internal 
           eig to compute the eigenvalues of B_k.           

Range of \hat\rho(k):
eigrange.m: A Matlab program that computes the range of \hat\rho(k).
            Usage: [minxi,maxxi] = eigrange(n,k);
            minxi is the minimum of \hat\rho(k) and maxxi is the maximum
            of \hat\rho(k).

Moments of \hat\rho(k):
mom.m: A Matlab program that computes the first five noncental and
       central moments of \hat\rho(k).
       The computation is based on explicit expressions of the moments
       of \hat\rho(k).
       Usage: [y1,y2] = mom(n,k);
       y1 is the first five noncentral moments of \hat\rho(k), y2 is the 
       first five central moments of \hat\rho(k).
moms.m: Symbolic version of mom.m (requires symbolic toolbox).
momrho.m: A Matlab program that computes the first s noncentral and
          central moments of \hat\rho(k). The computation is based on 
          a recursive formula.
          Usage: [y1,y2] = momrho(n,k,s);
          y1 is the first s noncentral moments of \hat\rho(k), y2 is the 
          first s central moments of \hat\rho(k).
momlinrho.m: A Matlab program that computes the first s noncentral and         
             central moments of 
             X = w_0+w_1*\hat\rho(1)+...+w_m*\hat\rho(m).
             The computation is based on a recursive formula.  The 
             program calls eigpsp.m to compute the eigenvalues of 
             a symmetric Toeplitz matrix.
             Usage: [y1,y2] = momlinrho(n,w,s);
             w is a vector of weights, y1 is the first s noncentral 
             moments of X, y2 is the first s central moments of X.
       
Cross-Moments of \hat\rho(k):
covrho.m: A Matlab program that computes Cov[\hat\rho(k_1),\hat\rho(k_2)].
          Usage: y = covrho(n,k1,k2);
covrhos.m: Symbolic version of covrho.m
covrho2.m: A Matlab program that computes Cov[\hat\rho(k_1)^2,\hat\rho(k_2)^2].
           Usage: y = covrho2(n,k1,k2);
covrho2s.m: Symbolic version of covrho2.m
pmomrho.m: A Matlab program that computes 
           E[\hat\rho(k_1)^s_1... \hat\rho(k_p)^s_p].
           Usage: y = pmomrho(n,k,s);
           k is a vector of [k_1, k_2, ..., k_p], s is a vector of
           [s_1, s_2, ..., s_p].

PDF of \hat\rho(k):
rhopdfb.m: A Matlab program that computes the pdf of \hat\rho(k).
           The program calls eigpapb.m to obtain the eigenevalues
           of B_k and then compute the pdf using a simplified version
           of the algorithm of Broda and Paolella (2009).
           Usage: f = rhopdfb(c,n,k);
           f is the density of \hat\rho(k) at c;
           For example, rhopdfb(0.4,20,4) returns 0.1444689987211.

CDF of \hat\rho(k):
rhocdfb.m: A Matlab program that computes the cdf of \hat\rho(k).
           The program calls eigpapb.m to obtain the eigenevalues
           of B_k and then linchi2b.m to compute the cdf of
           a linear combination of chi-squared random variables. 
           Usage: p = rhocdfb(c,n,k);
           p is P[\hat\rho(k)<c];
           For example, rhocdfb(0.4,20,4) returns 0.990830373693520.
irhocdfb.m: A Matlab program that computes the inverse cdf of \hat\rho(k).
            The program calls eigpapb.m to obtain the eigenevalues
            of B_k and then linchi2b.m to compute the cdf of
            a linear combination of chi-squared random variables. 
            Usage: c = irhocdfb(p,n,k);
            c is the constant such that P[\hat\rho(k)<c] = p.
            For example, irhocdfb(0.95,20,4) returns 0.273130088634109.

PDF of a linear combination of \hat\rho(k):
linrhopdf.m: A Matlab program that computes the pdf of 
             X = w_0+w_1*\hat\rho(1)+...+w_m*\hat\rho(m).
             The program calls eigpsp.m to compute the eigenvalues of
             a symmetric matrix.
             Usage: p = linrhopdf(c,n,w);
             w is a vector of weights, p is f_X(c).

CDF of a linear combination of \hat\rho(k):
linrhocdf.m: A Matlab program that computes the cdf of 
             X = w_0+w_1*\hat\rho(1)+...+w_m*\hat\rho(m).
             The program calls eigpsp.m to compute the eigenvalues of
             a symmetric matrix.
             Usage: p = linrhocdf(c,n,w);
             w is a vector of weights, p is P[X<c].
ilinrhocdf.m: A Matlab program that computes the inverse cdf of 
              X = w_0+w_1*\hat\rho(1)+...+w_m*\hat\rho(m).
              The program calls eigpsp.m to compute the eigenvalues of
              a symmetric matrix.
              Usage: c = ilinrhocdf(p,n,w);
              w is a vector of weights, c is the constant such that
              p = P[X<c].

Exact mean and variance of various autocorrelation based tests:
knoke.m: A Matlab program that computes the mean and variance of Knoke's (1977)
         statistic. 
         Usage: [y1,y2] = knoke(n);
         y1 and y2 are the mean and variance of Knoke's statistic,
         respectively.
vratio.m: A Matlab program that computes the mean and variance of the variance
          ratio test statistic \hat\theta(m).
          Usage: [y1,y2] = vratio(n,m);
          y1 and y2 are the mean and variance of \hat\theta(m), respectively.
lhr.m: A Matlab program that computes the mean and variance of the long-horizon
       regression test statistic \hat\beta(m).
       Usage: [y1,y2] = lhr(n,m);
       y1 and y2 are the mean and variance of \hat\beta(m), respectively.
bp.m: A Matlab program that computes the mean and variance of the Box-Pierce
      test statistic Q_{BP}(m).
      Usage: [y1,y2] = bp(n,m);
      y1 and y2 are the mean and variance of Q_{BP}(m), respectively.
lb.m: A Matlab program that computes the mean and variance of the Ljung-Box
      test statistic Q_{LB}(m).
      Usage: [y1,y2] = lb(n,m);
      y1 and y2 are the mean and variance of Q_{LB}(m), respectively.

Unadjusted and adjusted Box-Pierce test statistic:
qbp.m: A Matlab program that computes Q_{BP}(m) and Q_{BP}^a(m).
       Usage: [Q1,Q2] = qbp(x,m);
       x is a vector of the time series data, m is the number of lags,
       Q1 is the unadjusted Box-Pierce test statistic, Q2 is the
       adjusted Box-Pierce test statistic, with mean m and variance 2m.
