% Implementation of AJ Spectrogram: Plotting Histograms
% for SJ, SFA, SW, QVSplit and Beta
% Ait-Sahalia, Y. and Jacod, J., "Analyzing the Spectrum of Asset Returns:
% "Jump and Volatility Components in High Frequency Data"

% Characteristics of the data series
% Values below match the sample files for INTC and MSFT 2006

% read the  values of B(p,u,delta) computed and saved by the code
% AJspectrogram.m

year = '06';

stocknamevec = ['INTC','MSFT'];
qtrvec = [1;2;3;4];
Nstock = length(stocknamevec)/4; Nqtr = length(qtrvec);

SJ = 0; SFA = 0; SW = 0; QVSplit = 0; Beta = 0; % Initialize vectors of results


for stockindex=1:Nstock; % loop on stocks
    stockname = stocknamevec((1+(stockindex-1)*4):(stockindex*4));
    % read the stock names in successive blocks of 4 characters;
    
    for qtrindex=1:Nqtr; % loop on values of qtr
        qtr = qtrvec(qtrindex);
        
        % read saved values of B(p,u,delta)
        datafilename = [stockname,'_CleanTrans_5secFixed_',year,'Q',num2str(qtr),'_Results_B.mat'];
        load(datafilename,'-mat');
        % that mat file contains the variables B(pindex,aindex,gindex,dindex,kindex)
        % as well as pvec,atruncvec,gammavec,deltavec,kvec
        Np = length(pvec); Natrunc = length(atruncvec); Ngamma = length(gammavec);
        Ndelta = length(deltavec); Nk = length(kvec);
        
        if stockindex==1 && qtrindex==1; %initialize arrays based on length of pvec, etc., read in the first file
            U = zeros(Nstock,Nqtr,Natrunc,Ngamma,Ndelta);
        end;
        
        % assign in array the U's for use in Beta at p=0 only
        % recall that U = B(without truncation) - B(with truncation)
        for aindex=1:length(atruncvec);
            for gindex=1:length(gammavec);
                for dindex=1:length(deltavec);
                    U(stockindex,qtrindex,aindex,gindex,dindex) = B(1,Natrunc,gindex,dindex,1) - B(1,aindex,gindex,dindex,1);
                end;
            end;
        end;
        
        % compute SJ
        pvec_SJ = find((pvec(:)>=2.5 & pvec(:)<=6) );
        deltavec_SJ = find( deltavec(:)<=120 );
        kvec_SJ = find( kvec(:)>=2 );
        SJtemp=zeros(length(pvec_SJ),length(deltavec_SJ),length(kvec_SJ));
        for pindex=1:length(pvec_SJ);
            for dindex=1:length(deltavec_SJ);
                for kindex=1:length(kvec_SJ);
                    SJtemp(pindex,dindex,kindex) = B(pvec_SJ(pindex),Natrunc,1,deltavec_SJ(dindex),kvec_SJ(kindex)) / B(pvec_SJ(pindex),Natrunc,1,deltavec_SJ(dindex),1);
                end;
            end;
        end;
        SJtemp = reshape(SJtemp,prod(size(SJtemp)),1); % make the array into a column vector
        SJ = [SJ;SJtemp]; % append the new results to those for the other stocks/qtr
        
        
        % compute SFA
        pvec_SFA = find( (pvec(:)>=2.5 & pvec(:)<=6) );
        atruncvec_SFA = find( (atruncvec(:)>5 & atruncvec(:)<10) );
        deltavec_SFA = find( deltavec(:)<=120 );
        kvec_SFA = find( kvec(:)>=2 );
        SFAtemp=zeros(length(pvec_SFA),length(atruncvec_SFA),length(deltavec_SFA),length(kvec_SFA));
        for pindex=1:length(pvec_SFA);
            for aindex=1:length(atruncvec_SFA);
                for dindex=1:length(deltavec_SFA);
                    for kindex=1:length(kvec_SFA);
                        SFAtemp(pindex,aindex,dindex,kindex) = B(pvec_SFA(pindex),atruncvec_SFA(aindex),1,deltavec_SFA(dindex),kvec_SFA(kindex)) / B(pvec_SFA(pindex),atruncvec_SFA(aindex),1,deltavec_SFA(dindex),1);
                    end;
                end;
            end;
        end;
        SFAtemp = reshape(SFAtemp,prod(size(SFAtemp)),1); % make the array into a column vector
        SFA = [SFA;SFAtemp]; % append the new results to those for the other stocks/qtr
        
        
        % compute SW
        pvec_SW = find( (pvec(:)<2 & pvec(:)>1) );
        atruncvec_SW = find( (atruncvec(:)>5 & atruncvec(:)<10) );
        deltavec_SW = find( deltavec(:)<=120 );
        kvec_SW = find( kvec(:)>=2 );
        SWtemp=zeros(length(pvec_SW),length(atruncvec_SW),length(deltavec_SW),length(kvec_SW));
        for pindex=1:length(pvec_SW);
            for aindex=1:length(atruncvec_SW);
                for dindex=1:length(deltavec_SW);
                    for kindex=1:length(kvec_SW);
                        SWtemp(pindex,aindex,dindex,kindex) = B(pvec_SW(pindex),atruncvec_SW(aindex),1,deltavec_SW(dindex),1) / B(pvec_SW(pindex),atruncvec_SW(aindex),1,deltavec_SW(dindex),kvec_SW(kindex));
                    end;
                end;
            end;
        end;
        SWtemp = reshape(SWtemp,prod(size(SWtemp)),1); % make the array into a column vector
        SW = [SW;SWtemp]; % append the new results to those for the other stocks/qtr
        
        
        % compute QVSplit
        pvec_QVSplit = find( pvec(:)==2 );
        atruncvec_QVSplit = find( (atruncvec(:)>2 & atruncvec(:)<5) );
        deltavec_QVSplit = find( deltavec(:)<=1800 );
        QVSplittemp=zeros(length(atruncvec_QVSplit),length(deltavec_QVSplit));
        for aindex=1:length(atruncvec_QVSplit);
            for dindex=1:length(deltavec_QVSplit);
                QVSplittemp(aindex,dindex) = B(pvec_QVSplit(1),atruncvec_QVSplit(aindex),1,deltavec_QVSplit(dindex),1) / B(pvec_QVSplit(1),Natrunc,1,deltavec_QVSplit(dindex),1);
            end;
        end;
        QVSplittemp = reshape(QVSplittemp,prod(size(QVSplittemp)),1); % make the array into a column vector
        QVSplit = [QVSplit;QVSplittemp]; % append the new results to those for the other stocks/qtr
        
        
    end %loop on all 4 values of qtr
    
    
    % compute Beta
    % for Beta, aggregate all the year's data, otherwise not enough in one quarter
    if stockindex==1; %initialize Usum the first time around only
        Usum=zeros(Nstock,length(atruncvec),length(gammavec),length(deltavec));
    end;
    
    for aindex=1:length(atruncvec);
        for gindex=1:length(gammavec);
            for dindex=1:length(deltavec);
                % sum the number of large increments over the 4 quarters
                Usum(stockindex,aindex,gindex,dindex) = sum(U(stockindex,:,aindex,gindex,dindex));
            end;
        end;
    end;
    
    pvec_Beta = find( pvec(:)==0 );
    atruncvec_Beta = find( (atruncvec(:)>5 & atruncvec(:)<8) );
    gammavec_Beta = find( gammavec(:)>=1.5 & gammavec(:)<2.5 );
    deltavec_Beta = find( deltavec(:)<=10 );
    Betatemp=zeros(length(atruncvec_Beta),length(gammavec_Beta),length(deltavec_Beta));
    
    disp(length(atruncvec_Beta));
     disp(length(gammavec_Beta));
     disp(length(deltavec_Beta));
    
    for aindex=1:length(atruncvec_Beta);
        for gindex=1:length(gammavec_Beta);
            for dindex=1:length(deltavec_Beta);
                Betatemp(aindex,gindex,dindex) = (1/log(gammavec(gammavec_Beta(gindex)))) * log(Usum(stockindex,atruncvec_Beta(aindex),1,deltavec_Beta(dindex))/Usum(stockindex,atruncvec_Beta(aindex),gammavec_Beta(gindex),deltavec_Beta(dindex)));
            end;
        end;
    end;
    
    Betatemp = reshape(Betatemp,prod(size(Betatemp)),1); % make the array into a column vector
    Beta = [Beta;Betatemp]; % append the new results to those for the other stocks/qtr
    
    
    
end % loop on stock name (INTC and MSFT)






SJ = SJ(SJ(:)>0); SFA = SFA(SFA(:)>0); SW = SW(SW(:)>0);
QVSplit = QVSplit(QVSplit(:)>0); Beta = Beta(Beta(:)>0);
% eliminate the spurious 0 used to initialize in first place
% now display histograms

subplot(2,3,1), hist(SJ,0:0.1:3), title('SJ'), xlim([0 3]);
subplot(2,3,2), hist(SFA,0:0.1:3), title('SFA'), xlim([0 3]);
subplot(2,3,3), hist(SW,0:0.1:3), title('SW'), xlim([0 3]);
subplot(2,3,4), hist(QVSplit,0:0.1:1), title('QVSplit'), xlim([0 1]);
subplot(2,3,5), hist(Beta,0:0.2:2), title('\beta'), xlim([0 2]);
subplot(2,3,6), text(0.05,0.9,'SJ: test for jumps'),
text(0.15,0.8,'limit=1 if jumps, k^{p/2-1} if no jumps'),
text(0.05,0.7,'SFA: test for finite jump activity'),
text(0.15,0.6,'limit=1 if infinite activity, k^{p/2-1} if finite'),
text(0.05,0.5,'SW: test for presence of Brownian motion'),
text(0.15,0.4,'limit=1 if no Brownian, k^{1-p/2} if Brownian'),
text(0.05,0.3,'QVSplit: % of QV due to Brownian'),
text(0.05,0.2,'\beta: Estimate of degree of jump activity'), axis off ;








