% Implementation of AJ Spectrogram: Computation of B(p,u,delta)
% Ait-Sahalia, Y. and Jacod, J., "Analyzing the Spectrum of Asset Returns:
% "Jump and Volatility Components in High Frequency Data"

% this code computes the values of B(p,u,delta) and saves them
% run it only once, then do the plotting using AJspectrogram_Plots.m

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


% first read the data in the ascii file
% file format: DAY        TIME      TIME INTERVAL     PRICE     LOG-RETURN
% example:  03JAN2006	9:30:10          5            25.21     0.000396747

year = '06';

stocknamevec = ['INTC','MSFT'];

qtrvec = [1;2;3;4];

% will loop to compute the values of B(p,u,delta) for all those values below
pvec = [0;0.25;0.5;0.75;1;1.25;1.5;1.75;2;2.25;2.5;2.75;3;3.25;3.5;3.75;4;4.25;4.5;4.75;5;5.25;5.5;5.75;6];
atruncvec = [2;3;4;5;6;7;8;9;10;11;12;13;14;15;16;17;18;19;20;25;30;40;50;60;75;100;10^10];
% atrunc is expressed in terms of numbers of stdev of the continuous part
% use 10^10 to mean no truncation
gammavec = [1;1.25;1.5;1.75;2;2.25;2.5;2.75;3];

deltavec = [5;10;15;30;45;60;120;300;600;1800];
% for simplicity of coding, make sure those are multiples of timeinterval = 5 in the dataset
kvec = [1;2;3];

Nstock = length(stocknamevec)/4; Nqtr = length(qtrvec);
Np = length(pvec); Natrunc = length(atruncvec); Ngamma = length(gammavec);
Ndelta = length(deltavec); Nk = length(kvec);
B = zeros(Np,Natrunc,Ngamma,Ndelta,Nk);



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 ascii datafile
        % datafilename = INTC_CleanTrans_5secFixed_06Q1.txt
        datafilename = [stockname,'_CleanTrans_5secFixed_',year,'Q',num2str(qtr),'.txt'];
        
        [day,time,timeinterval,price,logreturn] = textread(datafilename,'%s %s %f %f %f','headerlines',1);
        
        dX = logreturn; % both X and dX have length n
        x0 = log(price(1)); % initial value
        
        n = length(dX); % this is the sample size
        T = 1/4; %l ength of time in years -- this is fixed, the sample size n is T/delta
        % 1/252 for one day, 21/252 for one month, 1/4 for one quarter, 1 for one year
        
        for dindex=1:Ndelta; % loop on values of sampling interval
            
            disp(['dindex = ',num2str(dindex)]);
            
            nblagj = deltavec(dindex)/deltavec(1); % measured in seconds, need to divide by deltavec(1) because nblagj is the number of lags
            deltaj = deltavec(dindex)/(6.5*60*60*252); % measured in years, this is the value of delta corresponding to that jindex
            % nj = n / nblagj; % this is the sample size of the vector Xobsj
            
            % To generate a series of numbers from 10 to 50, incrementing by 5, use A = 10:5:50
            % A = 10    15    20    25    30    35    40    45    50
            
            % observed increments at that frequency            
            X = x0 + cumsum(dX); % do this instead of X=log(price) to avoid including the large overnight returns
            dXobsj = X((nblagj+1):nblagj:n) - X(1:nblagj:(n-nblagj)); % length(dXobsj) is equal to nj-1
            disp(['dindex = ',num2str(dindex),' nb obs = ',num2str(length(dXobsj))]);
            
            % adapt==1
            % use the small increments in dXobsj to estimate sigma (= volatility of the continuous part of the semimartingale)
            sigmahat = sqrt( (1/T) * sum( (abs(dXobsj).^2) .* ( abs(dXobsj) <= 3 * 0.30 * deltaj^(1/2) )) );
            % use the increments that are smaller than the finite cutoff 3 * 0.30
            % the factor (1/T) is there in the variance to annualize the estimator sigmahat
            disp(['sigmahat = ',num2str(sigmahat)]);
            
            for kindex=1:Nk; % loop on values of k
                k = kvec(kindex);
                nblagjk = nblagj * k;
                dXobsjk = X((nblagjk+1):nblagjk:n) - X(1:nblagjk:(n-nblagjk));
                % for all values of kindex, do the truncation based on the value deltaj
                % so that we are using the same cutoffs u_n for delta and 2*delta
                
                for aindex=1:Natrunc; % loop on values of atrunc set as number of standard deviations
                    atrunc = atruncvec(aindex);
                    % the cutoff will be  u_n = atrunc * sigmahat * deltaj^(1/2)
                    % yes, it's based on deltaj, not deltajk
                    
                    for gindex=1:Ngamma;  % loop on values of gamma set as multiple of atrunc
                        g = gammavec(gindex);
                        
                        for pindex=1:Np; % loop on values of the exponent, to construct the B's
                            p = pvec(pindex);
                            B(pindex,aindex,gindex,dindex,kindex)  = sum( (abs(dXobsjk).^p) .* ( abs(dXobsjk) < g * atrunc * sigmahat * deltaj^(1/2) ));
                        end % of loop on pindex
                             
                    end % of loop on gindex                 
                end % of loop on aindex
            end % of loop on kindex
        end % of loop on dindex (delta)
        
                
        loopstotal = Nstock*Nqtr; loopsdone = qtrindex + (stockindex - 1)*Nqtr;
        disp(['data file # ',num2str(loopsdone),' out of ', num2str(loopstotal)]);
        
       
  
        % save data in matlab .mat files to be loaded later in Mathematica
        % do not include datapath so it saves in the same folder as AJspectrogram.m
        OUTFILE = [stockname,'_CleanTrans_5secFixed_',year,'Q',num2str(qtr),'_Results_B.mat'];
        
        
        % Save variables for MATLAB 5 & 6 compatibility in all cases, otherwise can't read the new (post version 7) mat files in Mathematica
        % this requires adding the option '-V6' (at the end of the save statement) when using Matlab version 7 or higher,
        % but this option is not recognized in earlier Matlab versions, so can't use it hence the else part
        
        if str2double(strtok(version,'.')) >= 7
            save(OUTFILE,'n','T','pvec','atruncvec','gammavec','deltavec','kvec', 'B','-V6');
        else
            save(OUTFILE,'n','T','pvec','atruncvec','gammavec','deltavec','kvec', 'B');
        end
        
        disp([' ']); disp(['File    ',num2str(OUTFILE),'   saved.']);
        
    end %loop on all 4 values of qtr
    
end % loop on stock name (INTC and MSFT)

