These are programs for my paper "Exact Variance Ratio Test with
Overlapping Data"
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 mrank.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: 1/10/2011, initial release
Version 1.1: 2/20/2011
Add optk_ar1r.m and optk_ar1rwp.m
Change initial estimate for ivr_cdf.m, ivr_ar1r_cdf.m, and ivr_ar1rwp_cdf.m
Change the names of some of the functions
Version 1.2: 11/22/2011
Add vr_pdf.m, vr_ar1r_pdf.m, vr_ar1rwp_pdf.m
Update all the cdf programs to work with matrix input


vr.m: Sample variance ratio and t-statistic using the Lo and MacKinlay (1988) method 
      (this function can operate on a TxN matrix of returns)
      Usage: [vr,trat] = vr(r,k);
      r is a TxN matrix of returns
      k is the length of multi-period return
      vr is a 1xN vector of \hat\theta(k), one for each column of data in r
      trat is a 1xN vector of t-ratio for testing H_0: \theta(k)=0 using the Lo and 
      MacKinlay (1988) method.
varvr.m: Exact variance of sample variance ratio
         Usage: y = varvr(T,k);
         It returns Var[\hat\theta(k)] under the assumption that returns are i.i.d. normal.
covvr.m: Exact covariance between two sample variance ratios
         Usage: y = covvr(T,k1,k2);
         It returns Cov[\hat\theta(k1),\hat\theta(k2)] under the assumption that
         returns are i.i.d. normal.
skewkurt.m: Computes coefficient of skewness and coefficient of excess kurtosis 
            of a sample variance ratio 
            Usage: [y1,y2] = skewkurt(T,k);
            y1 is the coefficient of skewness for \hat\theta(k), y2 is the coefficient
            of excess kurtosis for \hat\theta(k), both are obtained under the 
            assumption that returns are i.i.d. normal.
skewkurts.m: Computes coefficient of skewness and coefficient of excess kurtosis 
             of a sample variance ratio (requires the symbolic toolbox)
momvr.m: Computes arbitrary number of noncentral and central moments of
         a sample variance ratio.
         Usage: [y1,y2] = momvr(T,k,s);
         y1 is the first s noncentral moments of \hat\theta(k), y2 is the 
         first s central moments of \hat\theta(k).
pmomvr.m: Computes E[\hat\theta(k_1)^{s_1}...\hat\theta(k_p)^{s_p}].
          Usage: y = pmomvr(T,k,s);
          k is a vector of [k_1, k_2, ..., k_p], s is a vector of [s_1, s_2, ..., s_p].
          The computation is based on a recursive formula from Kan (2008).   

Stock price follows a random walk (returns are normal i.i.d.)
vr_pdf.m: exact pdf of sample variance ratio.
          Usage: y = vr_pdf(vr,T,k);  
          Example: vr_pdf(1.2,60,3) = 1.1051046616 tells us that when T=60, k=3, 
          f_{\hat\theta(k)}(1.2) = 1.1051046616
vr_cdf.m: exact cdf of sample variance ratio.
          Usage: y = vr_cdf(vr,T,k);  
          Example: vr_cdf(1.2,60,3) = 0.843274373 tells us that when T=60, k=3, 
          P[\hat\theta(k)<1.2] = 0.843274373
ivr_cdf.m: exact inverse cdf of sample variance ratio.
           Usage: y = ivr_cdf(prob,T,k); 
           Example: ivr_cdf(0.95,60,3) = 1.341079657 tells us that when T=60, k=3, 
           P[\hat\theta(k)<1.341079657] = 0.95

Stock return follows an AR(1)            
vr_ar1r_pdf.m: exact pdf of sample variance ratio. 
               Usage: y = vr_ar1r_pdf(vr,T,k,phi); 
               Example: vr_ar1r_pdf(1.2,60,3,0.1) = 1.711681234525 tells us that when 
               T=60, k=3, phi=0.1, f_{\hat\theta(k)}(1.2) = 1.711681234525
vr_ar1r_cdf.m: exact cdf of sample variance ratio. 
               Usage: y = vr_ar1r_cdf(vr,T,k,phi); 
               Example: vr_ar1r_cdf(1.2,60,3,0.1) = 0.6310491098 tells us that when 
               T=60, k=3, phi=0.1, P[\hat\theta(k)<1.2] = 0.6310491098
ivr_ar1r_cdf.m: exact inverse cdf of sample variance ratio.
                Usage: y = ivr_ar1r_cdf(prob,T,k,phi); 
                Example: ivr_ar1r_cdf(0.95,60,3,0.1) = 1.49778537618 tells us that when 
                T=60, k=3, phi=0.1, P[\hat\theta(k)<1.49778537618] = 0.95

Stock price is the sum of an AR(1) component and a random walk component
phi: AR(1) coefficient
r: Var[RW]/Var[AR1], ratio of the variance of the innovations in random walk component
   over the one for the AR(1) component (setting r=0 means the stock price follows an
   AR(1) process)
vr_ar1rwp_pdf.m: exact pdf of sample variance ratio. 
                 Usage: y = vr_ar1rwp_pdf(vr,T,k,phi,r); 
                 Example: vr_ar1rwp_pdf(1.2,60,3,0.95,0.5) = 1.0073472716384 tells us that 
                 when T=60, k=3, phi=0.95, r=0.5, f_{\hat\theta(k)}(1.2) = 1.0073472716384
vr_ar1rwp_cdf.m: exact cdf of sample variance ratio. 
                 Usage: y = vr_ar1rwp_cdf(vr,T,k,phi,r); 
                 Example: vr_ar1rwp_cdf(1.2,60,3,0.95,0.5) = 0.8660178228 tells us that 
                 when T=60, k=3, phi=0.95, r=0.5, P[\hat\theta(k)<1.2] = 0.8660178228
ivr_ar1rwp_cdf.m: exact inverse cdf of sample variance ratio.
                  Usage: y = ivr_ar1rwp_cdf(prob,T,k,phi); 
                  Example: ivr_ar1rwp_cdf(0.95,60,3,0.95,0.5) = 1.31705119049 tells us that 
                  when T=60, k=3, phi=0.95, r=0.5, P[\hat\theta(k)<1.31705119049] = 0.95

Optimal k for varaince ratio test
optk_ar1r.m: Finds out the optimal k for variance ratio test when the stock return
             follows an AR(1) process.  The answer is always k^*=2.
optk_ar1rwp.m: Finds out the optimal k for (left-tailed) variance ratio test when 
               the stock price is the sum of an AR(1) component and a random walk 
               component.
               Usage [k,maxprob] = optk_ar1rwp(T,phi,r,kinit,p,cutoff);
               kinit: initial estimate of optimal k (optional)
               p: size of the test (default is 0.05) 
               cutoff: a vector of cutoff points for the test (optional)
               Example: [k,maxprob] = optk_ar1rwp(240,0.85,0.5,2,0.05) returns
               k = 18, maxprob = 0.50140442.

The following four functions are called by the programs above and you generally
do not need to call them directly:
ggeig.m: eigenvalues of the P'H'HP matrix 
pnull.m: pnull(n) creates an nx(n-1) orthonormal matrix with its columns orthogonal to 1_n 
pspar.m: computes the eigenvalues and eigenvectors of P'*S*P, where P is orthogonal to
         1_n, and S=c1*I_n-c2*A, where A is a KMS matrix
linchi2.m: exact cdf of a linear combination of indepedent \chi^2_1 random variables 
