
%% Estimation options

clc
clear

nlags=1000;             % Augmentation order, set =1000 for Schwert's rule
demeaning_Robinson1=1;  % Demeaning for volatility-based measures 0 (No demean),1 (consant),2 (trend),3 (constant+trend),4 (constant +quadratic trend)
demeaning_Robinson2=4;  % Demeaning for volume-based measures

%Select Range dv1xdv2 to run the joint test. 

%Serie 1
liminf1=0;
limsup1=1;

%Serie 2
liminf2=0;
limsup2=1;

%Step for grid search
stepi=0.01;

%Input Data
pathinput='F:\FicherosTrabajoMatlab\LongMemoryVolumeVolatility\DataDJ30_5mins_Sampling\JAE_Files\AAPL_Data'; %Change this to your data path

%Outpur Data
pathoutput='F:\FicherosTrabajoMatlab\LongMemoryVolumeVolatility\DataDJ30_5mins_Sampling\JAE_Files\'; %Change this to your output path


%% Loading data
[Data, B]=xlsread(pathinput);
stock=B{2,2};
Z1=Data(:,3);  %Volatility measure
Z2=log(Data(:,4));  %Log Trading measure
Serie1='logRangeGK';         
Serie2='logVLM';                       
  
 %% BRRT MULTIVARIATE JOINT TESTING 

d1v=liminf1:stepi:limsup1;
d2v=liminf2:stepi:limsup2; 
LM=zeros(length(d1v),length(d2v));

conta=0;
 for i=1:length(d1v)
     
     for j=1:length(d2v)
         d1=d1v(i);
         d2=d2v(j);
         Salida = BRRTtest(Z1,Z2,d1,d2,demeaning_Robinson1,demeaning_Robinson2,nlags);
         LM(i,j)=Salida.LM; 
         conta=conta+1;
         cadena=['Joint Testing. Progress: ' num2str(100*conta/(length(d1v)*length(d2v))) '%'];
         disp(cadena);
         
     end
 end
 
% Estimating d: smallest LM 

a=min(min(LM));
[b1, b2]=find(LM==a);

check=BRRTtest(Z1,Z2,d1v(b1),d2v(b2),demeaning_Robinson1,demeaning_Robinson2,nlags);

%% Enhancing Output Check
Salida=[d2v;LM]; Salida2=[NaN; d1v']; Salida=[Salida2 Salida]; %Stores the values of the LM test in the grid for inspection
check.Salida=LM;   

resusu=zeros(6,1);
 auxi=min(LM');
 pos=find(auxi<chi2inv(0.95,2));
 eval(['check.lower95CI' Serie1 '=d1v(pos(1));']);
 eval(['check.upper95CI' Serie1 '=d1v(pos(end));']);
 resusu(1)=d1v(pos(1));
 resusu(2)=d1v(pos(end));
 
 auxi=min(LM);
 pos=find(auxi<chi2inv(0.95,2));
 eval(['check.lower95CI' Serie2 '=d2v(pos(1));']);
 eval(['check.upper95CI' Serie2 '=d2v(pos(end));']);
 resusu(3)=d2v(pos(1));
 resusu(4)=d2v(pos(end));
 
 %Identifying region for which d1=d2 leads to LM<5.99 (i.e., common d
 
 commond=[];
 for i=1:length(d1v)
     pos=find(d1v(i)==d2v);
     if isempty(pos)==0
         j=find(d2v==d1v(pos));
         if LM(i,j)<chi2inv(0.95,2)
             commond=[commond; d1v(pos)];
         end
     end
             
  end
         
if isempty(commond)==1
    check.lower95CIcommond=[];
    check.upper95CIcommond=[];
else
    check.lower95CIcommond=commond(1);
    check.upper95CIcommond=commond(end);
    
end
if isempty(check.lower95CIcommond)==1
    resusu(5)=NaN;
else
    resusu(5)=check.lower95CIcommond;
end
if isempty(check.upper95CIcommond)==1
    resusu(6)=NaN;
else
    resusu(6)=check.upper95CIcommond;
end

%% Preparing summary output 

Resumen=cell(2,12);
value1=['dhat_' Serie1];
value2=['dhat_' Serie2];
value3=['lower95CI_' Serie1];
value4=['upper95CI_' Serie1];
value5=['lower95CI_' Serie2];
value6=['upper95CI_' Serie2];

Resumen(1,:)={'Stock',value1,value3, value4,value2,value5,value6,'lower95common','upper95common','Correlation','K-order augment.' , 'Observations'};
Resumen(2,1)=B(2,2);
resusu=resusu';
V=corrcov(check.Sigmae);
if nlags==1000
    nlags = floor(4*(length(Z1)/100)^0.25);
    SR=1;
end
resusu2=[check.d(1) resusu(1:2) check.d(2) resusu(3:end) V(1,2) nlags];
Resumen(2,2:11)=num2cell(resusu2);
Resumen(2,12)=num2cell(length(Z2));

cadenaresumen=['ResultsLongMemoryTesting_' stock '.xlsx'];
Cadena3=[pathoutput cadenaresumen];
xlswrite(Cadena3,Resumen);
clc;
disp(Resumen)
if SR==1
     disp(['Note: augmentation based on Schwert Rule with nlags= 4*(T/100)^0.25=' num2str(nlags)]);
 
end
disp('dhat shows point estimates of d given by the d values in the grid that minimize the LM test statistic.')
disp('Lower and upper 95 are the lower and upper bounds of the projection of the 95% confidence ellipsoides on the axes.')
disp('Lower and upper common show the values of a common d contained in the 95% confidence ellipsoide (if not empty.')
disp('Correlation is the endogeneous correlation of estimated innovations.')