varatio.m

function [VR,Zk,Zhk,p_Zk,p_Zhk]=varatio(x,k)
% 
% Syntax: [VR,Zk,Zhk,p_Zk,p_Zhk]=varatio(x,k)
%
% Calculates the Variance Ratio Test (VR Test) of a time series x, with
% and without the heteroskedasticity assumption.
% input  : x is a vector of time serie (observed prices)
%          k is a scalar 
% output : VR is the value of the VR Test
%          Zk is the homoscedastic statistic of the variance ratio
%          Zhk is the heteroscedastic statistic of the variance ratio
%          p_Zk is the probability of rejection of the null Hypothesis H0 
%               that the time series x is a random walk under the
%               homoscedastic assumption while it's true.
%          p_Zhk is the probability of rejection of the null Hypothesis H0 
%               that the time series x is a random walk under the
%               heteroscedastic assumption while it's true.
%
% i.e : if p-value<=0.05 then one can reject H0 at the 5% level 
% *****************************************************************
% Reference:
% Lo A, MacKinley AC (1989): "The size and power of the variance ratio
% test in finite samples". Journal of Econometrics 40:203-238.
%
% Elaborated by : BEN HASSEN Anis 
% "Institut Suprieur de Gestion de Tunis" (ISG Tunis)
% University of Tunis
% 41, rue de la Libert - Cit Bouchoucha - C.P. : 2000 Le Bardo
% Tunisia
% University e-mail: http://www.isg.rnu.tn/
% Personal e-mail: [e-mail]
% _________________________________________________________________
% January 14, 2013.



%%
rt1=diff(log(x)); % one period return series
T=length(rt1); % length of one period return series
moy=mean(rt1); % mean of one period return series
v=var(rt1); % variance of one period return series
%%
i=1:T-k+1;
M=zeros(T-k+1,k);
for j=1:k,
    M(:,j)=i+j-1;
end
rtk=sum(rt1(M'))'; % k periods rate of return series
%%
m=k*(T-k+1)*(1-k/T);
VR=1/m*sum((rtk-k*moy).^2)/v; % Variance ratio statistic
%%
Zk=sqrt(T)*(VR-1)*(2*(2*k-1)*(k-1)/(3*k))^(-.5); % homoscedastic statistic
%%
j=1:k-1;
vec1=(2/k*(k-j)).^2;
rst=(rt1-moy).^2;
aux=zeros(1,k-1);
for i=1:k-1,
   aux(i)=rst(i+1:T)'*rst(1:T-i);
end
vec2=aux/((T-1)*v)^2;
Zhk=(VR-1)*(vec1*vec2')^(-.5); % heteroscedastic statistic
%%
p_Zk = 1- diff(normcdf([-abs(Zk) abs(Zk)],0,1)); %  % Two-tailed normal test for Zk
p_Zhk = 1- diff(normcdf([-abs(Zhk) abs(Zhk)],0,1)); %  % Two-tailed normal test for Zhk

% result = {'Variance ratio','homostat','p-value','heterostat','p-value';VR,Zk,p_Zk,Zhk,p_Zhk}