clear
format long;
tic;

load residuals.txt; % Load VAR(4) residual time series here.

data = abs(residuals(:,1).^2); %ANDREOU AND GHYSELS (2002): CUSUM on modulus/squared returns, choose power 1 or 2


N = size(data,1);

K = 12; % Set minimal period length 
L = 4;  % Set maximum amount of change points

K = floor(min(N/2,K)); 
L = floor(min(L, N/K -1));


LSE = zeros(L + 1,1);
BIC = zeros(L + 1,1);
AIC = zeros(L + 1,1);

demeaneddata=data-mean(data);
LSE(1) = demeaneddata'*demeaneddata;
BIC(1) = -N/2*log(LSE(1)/N);
AIC(1) = -N/2*log(LSE(1)/N);


CP = zeros(L+1, L+1);



for l = 1:L
    
    SSS = zeros(l+1,1);
    SSE = zeros(l+1,1);
    for ll = 1:l+1
        SSS(ll) = 1 + (ll-1)*K;
        SSE(ll) = ll*K;     
    end
    SSE(l+1) = N;    
    
    maxSSS = zeros(l+1,1);
    for ll = 2:l+1
        maxSSS(ll) = N + 1 - (l + 2 - ll)*K;
    end
    
    
    for ll = 1:l
        CP(l+1,ll) = SSE(ll);
    end
    SE = 0;
        for ll = 1:l+1
            submean = mean(data(SSS(ll):SSE(ll)));
            demeaneddata = data(SSS(ll):SSE(ll))-submean;
            SE = SE + demeaneddata'*demeaneddata;
        end   
    LSE(l+1) = SE;
    BIC(l+1) = N*log(LSE(l+1)/N)+l*log(N);
    AIC(l+1) = N*log(LSE(l+1)/N)+2*l;

    stopCrit = 0;

    while stopCrit < 0.5
        
        % Compute sum of squared errors SE
        SE = 0;
        %[SSS, SSE]
        for ll = 1:l+1
            submean = mean(data(SSS(ll):SSE(ll)));
            demeaneddata = data(SSS(ll):SSE(ll))-submean;
            SE = SE + demeaneddata'*demeaneddata;
        end        
        
        % Check for new LSE estimate
        if SE < LSE(l+1)
            LSE(l+1) = SE;
            BIC(l+1) = -N/2*log(LSE(l+1)/N)-l*0.5*log(N);
            
            for ll=1:l
                CP(l+1,ll) = SSS(ll+1);
            end
        end
        
        % Set new subsample
        subCrit = 0;
        ll = l;
        while subCrit < 0.5
            SSE(ll) = SSE(ll)+1;
            SSS(ll+1) = SSS(ll+1) + 1;
            if SSS(ll+1) > maxSSS(ll+1)
                ll = ll - 1;
                if ll < 0.5
                    subCrit = 1;
                    stopCrit = 1;
                end
            else 
                for lll = ll+1:l
                    SSE(lll) = SSE(ll)+(lll-ll)*K;
                    SSS(lll+1) = SSS(ll+1) + (lll-ll)*K;
                end
                subCrit = 1;
            end
        end
        
    end
    
end

N
K
L

CP
BIC
LSE
AIC

toc;

x=0:L;
figure;
plot(x,BIC);
